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I.  INTRODUCTION  AND  SUMMARY 


1.1  INTRODUCTION 

The  objective  of  the  Systems,  Science  and  Software  research 
program  is  to  extend  our  present  understanding  of  the  excitation  of 
seismic  waves  by  underground  explosions  and  earthquakes  in  order  to 
improve  the  United  States'  ability  to  monitor  compliance  with 
treaties  limiting  underground  nuclear  explosions.  This  report  is 
for  the  six  month  period  between  1  October,  1979,  and  1  April,  1980, 
and  includes  five  distinct  technical  sections.  These  are  briefly 
summarized  in  the  following  subsections. 

1.2  SUMMARY  OF  SECTION  II:  "A  SIMULATION  STUDY  OF  the 

DETECTABILITY  OF  A  5.3  KT  DECOUPLED  EXPLOSION  AT  REGIONAL 
DISTANCES  IN  THE  EASTERN  UNITED  STATES" 

This  section  was  written  by  J.  R.  Murphy  and  T.  J.  Bennett  of 
the  Reston  Geophysics  office  of  Systems,  Science  and  Software.  The 
objective  is  to  extend  the  regional  seismic  data  base  for  explosions 
in  the  eastern  United  States  by  theoretically  scaling  observed 
SALMON  data  to  simulate  the  seismograms  expected  from  a 

(hypothetical)  fully  decoupled  5.3  KT  explosion  at  the  SALMON  site. 
The  vertical  SALMON  recordings  from  four  stations  (EUAL,  a  *  246  km; 
CPO,  a  »  623  km;  SLWV,  a  =.  1065  km;  WFMN,  a  =*  1437  km)  were 
convolved  with  an  operator  which  is  intended  to  represent  the  ratio 
of  the  source  spectrum  for  a  decoupled  5.3  KT  event  to  the  source 
spectrum  for  SALMON.  Analysis  of  the  resulting  seismograms 

indicates  that  such  a  decoupled  event  would  probably  be  detectable 
at  the  two  closer  stations,  but  not  at  the  two  stations  with  a  > 
1000  km.  At  more  quiet  stations,  the  Lg  phase  would  probably  be 

detected  at  all  four  stations,  but  the  initial  P  wave  would  be 
difficult  to  detect. 

1.3  SUMMARY  OF  SECTION  III:  "YIELD  ESTIMATION  FROM  SURFACE  WAVE 
RECORDINGS  OF  UNDERGROUND  EXPLOSIONS  —  A  REVIEW" 

In  December,  1979,  DARPA  solicited  contributions  from 
interested  scientists  to  a  state-of-the-art  assessment  of  seismic 
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yield  determination.  S-Cubed  scientists  submitted  nine  separate 
summaries: 

1.  Definitions  of  Body  Wave  Magnitude 

T.  C.  Bache 

2.  Definition  of  Surface  Wave  Magnitude 

T.  C.  Bache 

3.  Source  Theory  and  Observation  for  Surface  Waves 

T.  C.  Bache 

4.  Tectonic  Generation  of  Surface  Waves 

T.  C.  Bache 

5.  Effects  of  Source  Region  Properties  on  Surface 
Wave  Generation 

T.  C.  Bache 

6.  Effects  of  Attenuation  on  Surface  Waves 

T.  C.  Bache 

7.  Body  Wave  Coupling  Theory 

J.  R.  Murphy 

8.  Experimental  Data  on  Body  Wave  Coupling 

J.  R.  Murphy 

9.  Source  Coupling  for  Body  Waves  and  Surface  Waves 

J.  T.  Cherry 

In  most  cases  these  summaries  are  very  useful  as  concise 
reviews  of  past  work  that  has  appeared  in  many  reports.  Those 
numbered  3,  4,  5  and  6  are  concerned  with  various  aspects  of  the 
problem  of  estimating  yield  from  single  station  recordings  of 
surface  waves.  The  material  in  these  four  summaries  has  been 
condensed  still  further  into  a  single  review  which  appears  as 
Section  III  of  this  report. 


1.4  SUMMARY  OF  SECTION  IV:  "THE  CONTRIBUTION  OF  TWO-DIMENSIONAL 
SOURCE  EFFECTS  TO  THE  FAR-FIELD  SEISMIC  SIGNATURES  OF 
UNDERGROUND  NUCLEAR  EXPLOSIONS" 

In  this  section  we  reprint  the  "Abstract"  and  "Introduction 
and  Summary"  from  a  report  of  the  same  title  (Bache,  et  aU ,  1980) 
just  submitted  to  VSC/ARPA. 


1.5  SUMMARY  OF  SECTION  V:  "THREE-DIMENSIONAL  EARTHOUAKE  MODELING 

INCLUDING  NONLINEAR  RUPTURE  OYNAMICS" 

Presented  in  this  section  is  a  summary  of  the  earthauake 

modeling  work  done  by  Steven  M.  Day  and  presented  to  the  AFOSR 

Program  Review  in  May,  1980. 

1.6  SUMMARY  OF  SECTION  VI:  "AUTOMATED  MAGNITUDES,  %  AND  Ms" 

This  section  was  written  by  T.  C.  Bache  and  describes 

algorithms  for  computing  automated  spectral  body  and  surface  wave 
magnitudes.  These  magnitudes  are  a  byproduct  of  MARS  signal  pro¬ 
cessing,  which  has  previously  been  implemented  for  signal  detection 
and  earthauake/explosion  discrimination.  They  promise  to  give  more 
consistent  and  convenient  measures  of  signal  energy  than  currently 
used  time  domain  magnitudes. 

An  earlier  version  of  the  mb  algorithm  was  described  and 

tasted  as  an  indicator  of  seismic  yield  by  Bache  (1979).  The 
current  algorithm  essentially  follows  the  recommendations  of  Masso 
et  al_.  (1979)  and  Savino,  et  al_.  (1980),  including  corrections  to 
account  for  the  presence  of  seismic  noise.  The  Ms  algorithm  is 
much  the  same,  except  that  the  dispersion  character  of  the  surface 
wave  signal  of  interest  must  be  pre-specif ied  so  the  time  series  can 
be  processed  to  make  a  matching  arrival  appear  to  be  undispersed. 

The  mb  and  M$  algorithms  are  tested  by  apDlying  them  to 
seismograms  constructed  by  adding  various  levels  of  seismic  noise  to 
a  signal  of  known  magnitude.  Consistent  magnitude  estimates  are 
obtained,  even  when  the  signal  is  nearly  obscured  by  the 
superimposed  noise. 

Also  presented  in  Section  VI  are  the  mb  results  for  a  series 
of  RKON  recordings  of  Eurasian  earthquakes  and  explosions  from  the 
AI  data  set.  The  magnitude  algorithms  are  ready  for  application  to 
a  much  larger  data  base. 
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II.  A  SIMULATION  STUDY  OF  THE  DETECTABILITY  OF  A 
5.3  KT  OECOUPLEO  EXPLOSION  AT  REGIONAL  DISTANCES 
IN  EASTERN  UNITED  STATES 

We  Have  recently  completed  a  preliminary  analysis  in  which  we 
theoretically  scaled  several  SALMON  short-period  seismograms  ob¬ 
served  at  regional  distances  in  the  eastern  united  States  to  produce 
estimates  of  the  corresponding  seismograms  which  would  be  expected 
at  these  stations  from  a  fully  decoupled  5.3  KT  explosion  at  the 
SALMON  shot  point.  The  primary  purpose  of  this  investigation  has 
been  to  extend  the  existing  regional  seismic  data  base  for  explo¬ 
sions  in  the  eastern  United  States  to  support  on-going  theoretical 


investigations  of 

detection. 

regional 

phase 

generation,  propagation 

and 

Four  eastern 

United 

States 

stations  were  selected 

for 

analysis:  Eutaw,  Alabama  (EUAL,  a  »  246  km),  Cumberland  Plateau 
Observatory,  Tennessee  (CPO,  a  *  623  km).  Beck  ley.  West  Virginia 
(BLWV,  a  *  1065  km)  and  Wykoff,  Minnesota  (WFMN,  a  -  1437  km).  The 
locations  of  these  four  stations  are  indicated  by  sauares  on  Figure 
2.1,  which  provides  a  map  view  of  the  North  American  stations  whicn 
recorded  the  SALMON  event  (Jordan  et  al .,  1966).  The  observed 
vertical  component  SALMON  data  recorded  at  these  four  stations  have 
been  theoretically  scaled  (Murphy,  1966)  to  the  ground  motion  to  be 
expected  from  a  fully  decoupled  5.3  KT  explosion  at  the  SALMON 
shotpoint.  The  theoretical  decoupling  factor  (with  respect  to 
SALMON)  corresponding  to  a  simple  step  in  pressure  in  a  39  m  radius 
cavity  is  shown  in  Figure  2.2.  It  should  be  noted  that  the  ampli¬ 
tude  level  in  Figure  2.2.  has  been  arbitrarily  normalized  to  give  a 
maximum  decoupling  factor  (at  2.5  Hz)  of  100.  This  reflects  the 
fact  that  the  observations  from  both  the  COWBOY  and  STERLING  experi¬ 
ments  indicate  low  freouency  decoupling  efficiencies  which  are  sig¬ 
nificantly  lower  than  those  predicted  theoretically.  Thus,  the 
expected  magnitude  of  the  decoupling  is  still  not  well  defined  and 
the  levels  shown  in  Figure  2.2  should  be  regarded  as  being  uncertain 
within  about  a  factor  of  two.  The  freouency  dependence  of  the  de¬ 
coupling  factor,  on  the  other  hand,  depends  primarily  on  the  shape 
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Figure  2.1.  Map  locations  of  regional  SALMON  stations  (□)  used  in  the  decoupling  simulation 


of  the  SALMON  source  function  for  frequencies  below  about  10  Hz,  and 
this  seems  to  be  reasonably  well  constrained  by  both  the  free-field 
and  near-regional  seismic  data. 


The  observed  SALMON  vertical  component  seismograms  recorded  at 
the  four  selected  stations  are  shown  in  Figures  2.3  to  2.6  (top), 
where  they  are  compared  with  the  corresponding  simulated  decoupled 
seismograms  (bottom).  The  simulated  decoupled  seismograms  were  ob¬ 
tained  by  convolving  the  observed  SALMON  seismograms  with  an 
operator  defined  by  the  ratio  of  the  decoupled  to  SALMON  theoretical 
source  functions.  The  traces  shown  to  the  left  of  the  simulated  de¬ 
coupled  seismograms  correspond  to  samples  of  the  noise  recorded 
prior  to  SALMON  at  each  station,  amplified  so  that  they  can  be  com¬ 
pared  with  the  estimated  decoupled  signal  levels.  This  provides  a 
basis  for  assessing  the  detectability  of  the  hypothesized  decoupled 
explosions  at  these  stations  under  the  local  noise  conditions  pre¬ 
vailing  at  the  time  of  the  SALMON  experiment.  A  comparison  of  the 
upper  and  lower  traces  in  these  figures  reveals  that,  in  each  of 
these  cases,  the  character  of  the  simulated  decoupled  seismogram  is 
very  similar  to  that  of  the  corresponding  SALMON  seismogram.  This 
reflects  the  fact  that  the  decoupling  factor  is  essentially  constant 
over  the  narrow  passband  available  from  the  standard  short-period 
systems  used  to  record  SALMON.  It  can  be  concluded  from  Figures  2.3 
to  2.6  that,  relative  to  these  background  noise  conditions,  such  a 
decoupled  explosion  would  probably  be  detectable  at  stations  EUAL 
and  CPO  (assuming  some  improvement  in  the  siqnal-to-noise  ratio 
resulting  from  signal  conditioning),  but  not  at  BLWV  or  WFMN. 
However,  it  must  be  noted  that  the  recorded  noise  levels  prior  to 
SALMON  (i.e.,  peak-to-peak  values  at  1  Hz  of  about  100  mw  at  EUAL 
and  30  mw  at  BLWV  and  WFMN)  are  significantly  higher  than  those 
associated  with  a  modern,  quiet  site  (i.e.,  about  2  mu  peak-to-peak 
at  1  Hz).  With  respect  to  this  ideal  noise  background,  Lg  would 
probably  be  detectable  at  all  four  stations*,  but  detection  of  the 

*  Scaling  data  are  not  presently  available  for  CPO.  However,  the 
signal  level  can  be  roughly  approximated  assuming  nominal  Lg 
distance  attenuation  with  respect  to  the  signal  level  recorded  at 
EUAL. 


Figure  2.4.  Comparison  of  observed  SALMON  and  simulated  decoupled  vertical  component 
seismograms,  station  CPO. 


Figure  2,4.  Continued. 


Sainton 


Figure  2.5.  Continued. 


Figure  2.6.  Comparison  of  observed  SALMON  and  simulated  decoupled  vertical  component 
seismograms,  station  WFMN. 


Fiqure  2.6.  Continued 


Figure  2.6.  Continued 
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initial  P  wave  motion  would  be  Questionable  at  the  two  distant 
stations  (BLWV  and  WFMN)  using  the  available  narrowband  data. 
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The  above  results  apply  to  data  recorded  through  narrowband 
instruments  and,  therefore,  tell  us  little  about  freauency  dependent 
detectability.  We  have  tried  to  remedy  this  by  correcting  for  in¬ 
strument  response.  Our  ability  to  accomolish  this  is  limited  by  two 
factors;  (1)  the  sampling  rate  used  in  digitizing  the  analog  data 
and,  (2)  the  freauency  dependent  signal-to-noise  ratio.  The  avail¬ 
able  data  were  sampled  at  a  rate  of  20  samples/second,  which  limits 
us  to  frequencies  below  10  Hz.  Estimates  of  the  signal  and  noise 
spectra  for  the  four  stations  are  shown  in  Figures  2.7  to  2.14.  For 
each  station,  smoothed  signal  spectra  (solid)  corresponding  to  50 
second  P  and  Lg  windows  are  compared  with  the  smoothed  noise 
spectrum  (dashed)  computed  from  the  50  second  window  immediately 
preceding  the  first  arrival.  The  dotted  lines  on  these  figures 
correspond  to  estimates  of  the  decoupled  signal  spectra  obtained  by 
scaling  down  the  corresponding  SALMON  spectra,  using  the  decoupling 
factor  shown  in  Figure  2.2.  It  can  be  seen  that  the  SALMON  signal- 
to-noise  ratios  at  EUAL,  CPO,  and  BLWV  are  good  over  the  freauency 
range  from  0.5  to  10  Hz,  while  those  for  the  most  distant  station 
(WFMN)  fall  below  unity  at  both  high  and  low  frequencies.  Figures 
2.15  to  2.17  compare  the  uncorrected  and  selected  corrected  band- 
widths  for  the  four  stations,  assuming  a  nominal  short-period  in¬ 
strument  response  for  each.  (The  actual  system  responses  at  these 
stations  during  SALMON  varied  somewhat  and  therefore  the  derived 
instrument  correction  factors  are  only  approxiate).  The  instrument- 
corrected  ground  motions  for  the  four  stations  are  shown  in  Figures 
2.18  to  2.21  (center),  where  they  are  compared  with  the  original 
SALMON  recordings  (top)  and  the  simulated,  instrument-corrected 
ground  motions  from  the  5.3  KT  decoupled  explosion  (bottom).  Here 
the  accompanying  noise  samples  correspond  to  instrument-corrected 
versions  of  those  shown  in  Figures  2.3  to  2.6.  It  can  be  seen  that 
the  instrument  correction  significantly  changes  the  character  of  the 
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Relative  Spectral  Amplitude 


Frequency  (Hz) 

Figure  2.7.  Comparison  of  P  wave  and  noise  spectral 
amplitude  levels,  station  EUAL. 
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Frquency  (Hz) 


Figure  2.8.  Comparison  of  La  and  noise  spectral  amplitude 
levels,  stat1onyEUAl. 
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Figure  2.1 0.  Comparison  of  lg  and  noise  spectral  amplitude 
levels,  station  CPO. 
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Salmon 


Frequency  (Hz) 

Figure  2.14.  Comparison  of  lg  and  noise  spectral  amplitude 
levels,  station  WFMN. 
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Fiaure  2. IB.  Comparison  of  Instrument-corrected  SALMON  and  simulated  decoupled 
vertical  component  seismograms,  station  EUAL. 


Sainton 


Fiqure  2.19,  Comparison  of  instrument-corrected  SALMON  and  simulated  decoupled 
vertical  component  seismograms,  station  CPO. 


Fiqure  2.19.  Continued. 


Fiqure  2.20.  Comparison  of  Instrument-corrected  SALMON  and  simulated  decoupled 
vertical  component  seismograms,  station  BLWV. 
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Figure  2,20.  Continued 


Salmon 


Figure  2.21.  Comparison  of  Instrument-corrected  SALMON  and  simulated  decoupled 
vertical  component  seismograms,  station  WFMN. 
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Figure  2.21.  Continued. 
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recordings,  particularly  at  stations  EUAl  and  CPQ.  However,  com¬ 
parisons  with  Figures  2.3  to  2.6  indicate  that  the  effective  signal- 
to-noise  ratios  predicted  by  the  decoupling  simulation  actually 
decrease  after  correction  for  instrument  response.  This  reflects 
the  fact  that  the  noise  level  increases  with  decreasing  freauency 
while  the  decoupling  factor  remains  essentially  constant.  Thus,  in 
this  case,  it  is  the  high  freauency  information  which  is  of  greatest 
potential  utility  and  corrections  for  low  freauency  instrument 
response  characteristics  are  not  expected  to  improve  detection 
capability.  However,  because  of  the  limited  sampling  rate,  it  has 
not  been  possible  to  auantitatively  evaluate  high  freauency 
detectability  using  these  data. 

A  final  auestion  is  whether  any  implications  can  be  drawn  from 
these  data  concerning  either  the  SALMON  source  function  or  anelastic 
attenuation  in  the  eastern  United  States.  Figure  2.2?  shows  a 
comparison  of  the  instrument-corrected  (JM)  SALMON  P  wave  spectrum 
from  CPO  with  a  theoretical  P  wave  spectrum  obtained  by  multiplying 
an  analytic  approximation  to  the  SALMON  source  function  (Murphy, 
1969)  times  an  attenuation  operator,  e-wft*,  with  t*  «  0.075.  rt 
can  be  seen  that  the  agreement  between  the  two  is  auite  good  for 
freauencies  above  1  Hz.  Below  1  Hz  the  two  curves  diverge,  but  it 
seems  likely  that  this  is  due  to  the  relatively  long  time  window 
(i.,e.  50  seconds)  used  in  computing  the  P  wave  spectrum.  That  is, 
the  selected  window  includes  late-arriving,  lower  freauency 
components  which  attenuate  less  rapidly  than  the  first  arriving  p 
wave  phases.  The  duality  of  the  fit  above  1  Hz  suggests  both  that 
the  selected  SALMON  source  function  is  about  right,  and  that  the 
anelastic  attenuation  is  auite  low  for  this  particular  eastern 
United  States  travel  path.  In  fact,  with  t*  -  r/cQ,  where  r  is  the 
travel  path  length  and  c  the  velocity  of  the  phase  under  consider¬ 
ation,  it  follows  that  Q  *  r/ct*.  For  r  >  i  and,  at  this  range,  c  < 
8  km/sec,  and  it  follows  that  a  t*  value  of  0.075  corresponds  to  an 
average  path  0  value  of  at  least  1000. 

In  summary,  regional  SALMON  data  recorded  at  four  eastern 
United  States  stations  have  been  theoretical ly  scaled  to  simulate 
the  seismograms  to  be  expected  from  a  fully  decoupled  5.3  KT 
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Figure  2.22.  Comparison  of  theoretical  and  instrument-corrected 
observed  P  wave  spectra,  station  CPO. 


explosion  at  the  SALMON  shotpoint.  It  has  been  found  that,  relative 
to  the  noise  background  prevailing  during  SALMON,  such  an  explosion 
would  probably  have  been  detectable  on  the  short-period  recordings 
from  the  two  nearer  stations  (EUAL,  a  ■  246  km  and  CPO,  a  ■  623  km), 
but  not  from  those  expected  at  the  more  distant  stations  (BLWV,  a  « 
1065  km  and  WFMN,  a  -  1437  km).  With  respect  to  the  lower  noise 
background  expected  at  a  modern,  quiet  site,  the  Lg  phase  from  sucn 
a  decoupled  explosion  would  probably  be  detectable  at  all  four 
stations,  but  initial  P  wave  detection  would  still  be  questionable 
at  the  two  distant  stations.  Furthermore,  as  might  be  expected  on 
the  basis  of  previous  noise  studies,  increasing  the  effective 
bandwidth  of  the  short-period  data  to  encompass  lower  frequencies 
does  not  improve  this  detectability  oue  to  the  fact  that  the  noise 
level  increases  with  decreasing  frequency,  while  the  decoupling 
factor  remains  essentially  constant.  By  tne  same  argument, 
consideration  of  higher  frequency  data  should  improve  detectability, 
at  least  for  high  q  paths.  However,  it  has  not  been  possible  to 
quantitatively  evaluate  this  potential  improvement,  using  the 
available  data  set,  due  to  the  relatively  low  sampling  rate  employed 
in  the  digitization.  Therefore,  it  is  recommended  that  selected 
SALMON  data  be  redigitized  at  a  higher  rate  so  that  frequency 
dependent  detection  thresholds  can  be  assessed  for  these  eastern 
United  States  paths. 


39 


SYSTEMS,  SCIENCE  AND  SOFTWARE 


i* 


I 

THIS  PAGE  LEFT  BLANK 


III.  YIELD  ESTIMATION  FROM  SURFACE  WAVE  RECORDINGS 
OF  UNDERGROUND  EXPLOSIONS  —  A  REVIEW 


3.1  INTRODUCTION 

In  this  section  we  summarize  the  results  from  several  stuoies 
of  surface  waves  that  have  been  done  over  the  past  several  years. 

The  main  objective  in  this  work  has  been  to  understand  the  para¬ 

meters  that  control  single  station  recordings  of  the  surface  waves. 
This  is  said  to  contrast  our  approach  from  alternatives  in  which  the 
data  are  represented  by  network  quantities  like  M$. 

This  review  includes  four  topics  which  are  discussed  in  sep¬ 
arate  sub-sections.  First,  we  summarize  a  study  of  Airy  phase 

amplitudes  of  NTS  explosions  recorded  at  the  WWSSN  stations  ALQ  and 
TUC.  The  purpose  of  this  study  was  to  delineate  the  amplitude  oif- 
ferences  among  events  in  different  source  materials.  These  aiffer- 
ences  are  compatible  with  the  differences  seen  in  data,  where 

the  two  can  be  compared.  The  ALQ  and  TUC  observations  are  the 
primary  data  base  for  subsequent  synthetic  seismogram  studies. 

The  second  topic,  discussed  in  Section  3.3,  is  the  inference 
of  explosion  reduced  displacement  potential  (RDP)  amplitudes  from 
comparison  of  synthetic  and  observed  waveforms.  The  synthetics  are 
computed  with  an  RDP  source  in  plane-layereo  earth  moaels.  The 
capability  to  construct  synthetic  seismograms  with  such  moaels  is 
widely  available,  and  it  is  important  to  see  how  well  it  can  explain 
the  data.  Of  course,  the  connection  between  explosion  yielo  and  RDP 
amplitude  is  another  area  where  our  understanding  is  incomplete,  but 
this  is  not  a  topic  for  oiscussion  here. 

Attempts  to  use  theoretical  models  to  quantitatively  explain 
the  data  have  not  often  been  made.  The  study  of  ALQ  and  TUC  data 
summarized  in  Section  3.3  is,  perhaps,  the  most  complete.  This 
study  demonstrates  that  the  models  do  give  results  that  are  in 
fairly  good  agreement  with  the  data.  However,  the  comparison  of 
theoretical  and  observed  seismograms  also  highlights  those  features 
that  cannot  be  explained  with  simple  source  and  parr  models.  In 
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Section  3.4  we  discuss  the  most  important  higher  oroer  source  con¬ 
tributions  and  the  extent  to  which  their  effects  can  be  delineateo 
with  elastic  point  source  models.  While  the  seismogram  synthesis 
can  handle  equivalent  elastic  point  sources  of  arbitrary  complexity, 
the  problem  is  in  specifying  the  source  components.  To  do  so,  an 
improved  understanding  of  the  physics  of  explosions  in  two-  and 
three-dimensions  is  required.  The  capability  to  move  in  that 
direction  is  now  available,  as  evidenced  by  the  results  from 
two-dimensional,  axisymmetric  source  calculations  which  are 
discussed  in  Section  IV  of  this  report. 

The  final  topic  of  this  summary  is  discussed  in  Section  3.5. 
In  that  section  we  point  out  that  many  underground  explosions  are 
characterized  by  strong  radiation  patterns,  but  that  a  satisfactory 
explanation  for  the  origin  of  this  phenomenon  is  not  yet  available. 
The  conventional  idea  that  a  double-couple  source  is  superimposed  on 
the  explosion  qualitatively  explains  many  of  the  observations,  but 
quantitative  corrections  cannot  confidently  be  made  ana  many 
important  features  remain  unexplained.  These  points  are  illustrated 
by  a  summary  of  results  from  a  study  of  the  waves  from  MIGHTY  EPIC 
and  DIABLO  HAWK,  two  nearly  identical  events  at  Rainier  Mesa. 

We  conclude  this  review  with  a  summary  of  the  main  conclusions 
in  Section  3.6.  The  theoretical  techniques  are  availaole  to  model 
most,  if  not  all,  of  the  important  parameters  controlling  surface 
waves  and  there  are  ample  data.  More  systematic  modeling  studies 
are  needed  to  clearly  define  the  uncertainty  in  yield  estimation 
when  the  most  powerful  theoretical  techniques  are  employed. 

3.2  EXPERIMENTAL  DETERMINATION  OF  THE  DEPENDENCE  OF  RAYLEIGH  WAVE 

AMPLITUDE  ON  PROPERTIES  OF  THE  SOURCE  MATERIAL 

There  have  been  numerous  studies  in  which  Rayleigh  wave 
amplitude,  usually  indicated  by  M  ,  is  plotted  versus  explosion 
yield.  Events  are  often  separated  according  to  characteristics  of 
the  source  materials  and  some  dependence  is  sought,  we  have  done  a 
similar  study  using  the  Airy  phase  amplitudes  measured  on  recordings 
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at  the  kWSSN  stations  ALQ  and  TUC.  These  stations  are  at  ranges 
from  700  to  900  kilometers  and  record  small  events  not  usually 
present  in  data  sets.  As  a  disadvantage,  they  are  off-scale 
for  events  mucn  larger  than  300  KT. 

The  results  of  this  study  were  described  in  a  1977  Systems, 
Science  and  Software  (S'*)  Quarterly  Report  by  Bache,  Goupillauo 
and  Mason.  The  Airy  phase  data  were  compared  to  M  data  compiled 
by  Eisenhauer  (1976).  The  Airy  phase  amplitudes  were  converted  to 
M$  values  using  the  formulas: 

ALQ:  *  log  A  +■  2.72, 

TUC:  M^  a  log  A  +  2.17, 

where  the  constants  were  chosen  to  make  M^  and  about 

the  same,  on  the  average,  as  Eisenhauer’s  M  .  An  M  was  then 

AT  s  s 

taken  to  be  the  mean  of  M$  and  M$. 

The  data  were  separated  by  test  area  and  a  linear  least  square 
fit  was  made  for  each  population  when  plotted  versus  log  yield.  The 
linear  best-fit  equations,  together  with  the  standard  deviations, 
are  as  follows  (the  number  of  events  is  listed  in  parentheses  ana  B 
is  an  arbitrary  constant): 

Pahute  Mesa  below  the  water  taole  (9): 

M$  =•  Q.86  log  W  +  B,  a  *  0.09, 

Yucca  Flat  below  the  water  table  (30): 

M  «  1.17  log  W  +  8  -  0.93,  a  =*  0.19, 

Pahute  Mesa  above  the  water  table  (13): 

Ms  »  1.21  log  W  +  B  -  1.09,  a  *  0.05, 

Yucca  Flat  above  the  water  table  (15): 

M$  .  0.78  log  W  +  B  -  0.32,  a  *  0.27, 

Rainier  Mesa,  tunnel  shots  only  (10) 

Ms  .  0.77  log  W  +  B  -  0.41,  o  .  0.19. 


t 
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The  M$  data  together  with  the  best  fitting  lines  are  plotted 
in  Figure  3.1.  As  expected  from  the  large  a  values  for  some  popu¬ 
lations,  there  is  considerable  scatter  in  the  data. 

There  are  seven  Pahute  Mesa  and  eleven  Yucca  Flat  events  for 
which  the  Airy  phase  estimates  can  be  compared  to  the  tele- 

seismic  of  Eisenhauer  (1976).  These  are  plotted  in  Figure 

3.2.  To  give  some  idea  of  the  scale,  the  standard  deviation  of  the 
mean  residual  between  the  two  M$  measurements  is  less  than  0.10 
M$  units.  This  indicates  that  the  Airy  phase  measurements  from 
these  two  stations  give  Mg  measurements  nearly  the  same  as  the 
tele-  seismic  M?  from  many  stations  compiled  by  Eisenhauer. 

The  data  plotted  in  Figure  3.1  indicate  that  the  slope  of  M 
versus  log  W  is  not  much  different  from  unity.  Then  one  way  to  in¬ 
dicate  the  relative  coupling  in  different  areas  is  to  compute  M$  - 
log  W  for  each  event.  The  mean  values  of  this  quantity  are  shown  in 
Figure  3.3  for  each  population. 

We  see  that  the  M$  and  Airy  phase  data  agree  that  the  M$ 
coupling  is  about  0.1  -  0.3  units  higher  for  events  below  the  water 
table  at  Pahute  Mesa  than  for  comparable  events  at  Yucca  Flat. 

There  are  only  a  few  events  above  the  water  table  at  Pahute  Mesa, 

but  they  clearly  couple  more  weakly  into  M$  than  those  in 
saturated  materials.  The  lowest  M$  events  are  those  in  dry  tuffs 

at  Yucca  Flat.  The  scatter  is  ouite  large  for  these  low  yield 

events,  but  M$  -log  W  is  0.57  lower,  on  the  average,  than  for  the 
saturated  events  at  Yucca  Flat.  This  difference  is  much  larger  than 
the  standard  deviation  of  the  data.  The  saturated  tuff  explosions 
in  the  tunnel  beds  at  Rainier  Mesa  seem  to  couple  about  the  same  as 
the  saturated  tuff  events  identified  as  being  in  Paleozoic  rock  at 
Yucca  Flat.  Since  these  events  were  detonatea  close  to  the  tuff- 
Paleozoic  interface,  the  identification  of  these  events  as  being  in 
Paleozoic  rock  is  somewhat  ambiguous,  and  we  are  probably  not  seeing 
the  true  differences  between  events  in  the  two  source  media. 
Finally,  the  PILEDRIVER  event  seems  to  couple  like  the  highest 
coupling  population,  the  Pahute  Mesa  events  below  the  static  water 
table. 
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Figure  3.1.  The  M$  values  from  the  WWSSN  stations  ALQ  and  TUC  are  plotted  versus  explosion 
yield.  The  explosions  are  divided  Into  different  populations  according  to  the 
gross  properties  of  the  emplacement  media. 


represent  the  mean  +  one  standard  deviation. 


In  summary,  we  have  delineated  the  coupling  differences 
between  events  in  different  areas  at  NTS.  The  Airy  phase  amplitudes 

at  ALQ  and  TUC  show  about  the  same  variation  with  yield  and  source 

material  as  the  from  a  much  larger  data  base.  The  question  is, 
do  we  understand  why  the  surface  wave  amplitudes  vary  from  one  test 
area  to  another?  This  question  motivates  the  work  discussed  in  the 
next  several  sections. 

3.3  EXPLOSION  RDP  AMPLITUDES  FROM  RAYLEIGH  WAVES 

3.3.1  Introduction 

The  most  obvious  way  to  think  about  using  surface  waves  to 
determine  explosion  yields  is  in  terms  of  the  Ms-yield 

relationship.  But  is  an  average  of  values  from  many  single 

station  recordings,  so  a  close  look  at  the  single  station  behavior 
of  the  Rayleigh  wave  amplitude  should  enhance  our  understanding  of 
what  is  controlling  the  M  .  In  two  reports  (Bache,  Rodi  and 

Harkrider,  (1978),  and  Bache,  Rodi  and  Masso,  (1978),  we  looked 
closely  at  the  ALQ  and  TUC  recordings  of  NTS  explosions  that  we 
discussed  in  the  last  section. 

The  simplest  model  for  the  explosion  source  is  the  reduced 
displacement  potential  (ROP)  and  a  great  deal  of  work  has  been  done 
to  interpret  the  data  in  terms  of  this  spherically  symmetric  source 
representation.  If  the  RDP  is  viewed  as  a  first-order  model  for  the 
source,  plane-layered  structures  provide  comparable  first-order 
models  for  the  travel  path.  In  the  work  described  in  the  two  1978 
reports,  the  observed  seismograms  were  interpreted  using  these 
fairly  simple  models.  In  this  section  we  will  summarize  that  work. 

3.3.2  Path  Models  Compatible  with  Surface  Wave  Observations 

To  deduce  the  source  amplitude  from  far-field  seismic 
recordings,  it  is  necessary  to  correct  for  the  effect  of  the  path. 
Bache,  Rodi  and  Harkrider  (1973)  constructed  path  models  for  the 
NTS- ALQ  and  NTS-TUC  paths  from  observed  surface  waves  at  the  two 
stations.  THe  main  elements  of  that  study  were: 
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•  Rayleigh  wave  recordings  of  NTS  explosions  were 
collected  from  the  VkWSSN  stations  ALQ  and  TUC. 

The  data  were  divided  according  to  test  site.  It 
was  observed  that  recordings  from  a  single  test 
site  (e.g.,  Pahute  Mesa)  had  very  similar  wave¬ 
forms  . 

•  Representative  recordings  were  digitized  from 
three  sites;  Climax  Stock,  Yucca  Flat  and  Pahute 
Mesa. 

•j 

•  The  digitized  data  were  processed  by  the  S 
MARS  program  and  phase  and  group  velocity  curves 
were  derived  for  the  two  paths,  NTS-TUC  and  NTS- 
ALQ. 

•  Using  generalized  linear  inversion  techniaues, 
plane  layered  earth  models  were  found  that  fit 
the  data.  These  models  are  compatible  with  other 
geophysical  information  about  these  paths. 

•  A  Q  model  was  constructed  from  western  united 
States  attenuation  data  collected  by  Mitchell 
(1975). 

•  Synthetic  seismograms  were  computed  with  a  simple 
reduced  displacement  potential  (ROP)  source. 

These  seismograms  are  shown  compared  to  typical 
observations  in  Figure  3.4. 

This  procedure  ensures  that  the  dispersion  of  the  synthetic 
and  observed  seismograms  be  the  same.  However,  the  excellent  wave¬ 
form  agreement  indicates  that  the  amplitude  spectra  are  fit  as  well. 


49 


arsreMS  sc/r/vcr  and  sof-twa* r 


Figure  3.4.  Theoretical  and  observed  seismograms  are  compared 
at  ALO  (left)  and  TUC  for  events  In  three  test 
areas  at  NTS.  A  tar  Indicating  one  minute  is 
shown.  In  each  pair  the  observed  (top)  and  theo¬ 
retical  records  start  at  the  same  time  with  res¬ 
pect  to  the  explosion  detonation  and  this  time  is 
indicated  as  T  . 
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3.3.3  Surface  Wave  Dependence  on  Source  Material  Properties 

In  computing  the  seismograms  of  Figure  3.4  there  is  one  path 
model  for  NTS-TUC  and  one  for  NTS-ALO,  but  we  must  deal  with  the 
fact  that  the  local  material  properties  are  different  for  the  three 
test  areas.  This  was  handled  by  using  a  two  path  model  with  the 
local  differences  accounted  for  in  the  top  two  kilometers  of  the 
source  region  model.  Transition  between  the  two  is  accounted  for  by 
an  approximate  transmission  coefficient,  TU)»  based  on  results  of 
McGarr  (1969).  The  vertical  Rayleigh  wave  is  computed  from 


w(r,w)  »  -4iru, 


$U) 


Turk 


(2) 


(f2) 


,-vr 


1/2 


sin  a 


where  subscripts  1  and  2  indicate  the  source  and  path  models. 

An  important  question  is,  keeping  all  other  factors  fixed,  how 

does  Rayleigh  wave  amplitude  scale  with  the  source  material  proper- 

* 

ties?  The  Y'U)  is  the  RDP  source  and  is^at  long  period,  us  is 
the  shear  modulus  at  the  source,  is  the  depth  dependent  eigen¬ 
function  for  an  explosion,  c  is  phase  velocity  and  AR  is  depth-in¬ 
dependent  amplification. 

The  answer  to  the  question  is  shown  in  a  direct  way  in  Figure 
3.5.  For  source  regions  that  are  not  too  different  (e.g.,  Pahute 
Mesa  ana  Yucca  Flat)  we  have 


Ms  *  log  („s  *„). 

However,  if  we  compare  the  Climax  Stock  to  the  others,  we  see  that 
the  relationship  has  a  strong  frequency  dependence  and  takes  no 
simple  form. 


3.3.4  Inversion  for  the  RDP  Source  Amplitude 

We  have  pointed  out  that  plane-layered  earth  models  and  an  RDP 
source  are  adequate  to  give  synthetic  seismograms  with  waveforms 
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Figure  3.5. 


(c) 


(a)  The  density,  shear  and  congressional  wave  veloc¬ 
ities  are  plotted  versus  depth  for  the  three  test 
area  at  NTS.  (b)  The  source  amplification  factor 
is  shown  for  the  three  source  areas  studied,  (c)  The 
transmission  coefficient  T(u>)  is  plotted  for  the  six 
source-path  combinations  studied. 


that  match  the  observations.  Using  this  comparison,  we  inferred 
the  ’fao  reauired  to  match  the  observed  amplitude.  This  work  is 
described  in  an  S3  report  by  Bache,  Rodi  and  Mason  (1978) 

The  events  considered  were  separated  into  the  three  popu¬ 
lations  and  all  were  below  the  water  table  with  yields  between  40 
and  200  KT.  The  inferred  values  are  shown  in  Figure  3.6.  We 
conclude  that: 

•  For  each  station  the  (scaled  to  a  common  yield) 

was  consistent  within  each  population  -  the 
standard  deviation  was  about  40  percent  of  the 
mean. 

•  The  f,,  from  one  station  (TUC)  was  consistently  50 
percent  larger  than  from  the  other. 

•  The  values  were  consistent  with  those  derived 
using  other  methods  (close-in  observations,  tele- 
seismic  body  and  surface  waves,  finite  difference 
source  calculations). 

These  events  do  not  span  a  wide  yield  or  depth  range.  The 
only  depth  correction  made  was  for  the  depth  dependence  of  the 
eigenfunctions.  This  turns  out  to  be: 

Yucca  Flat:  log  A  «  0.16  log  H, 

Pahute  Mesa:  log  A  ~  0.05  log  H. 

If  H  a  W1'3,  and  the  source  function  scales  with  yield,  this 

translates  to  an  M$  -  log  yield  slope  of  1.05  for  Yucca  Flat  and 
1.02  for  Pahute  Mesa. 

We  have  not  put  forth  an  explanation  for  the  relative 

amplitudes  of  explosions  in  the  other  classes  listed  in  Figure  3.3. 
The  events  above  tne  water  table  have  generally  lower  u  ana  f^than 

those  below  the  water  table,  so  the  observed  coupling  of  these 
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events  is  as  expected.  The  Rainier  Mesa  events  have  v  close  to 
that  for  Yucca  Flat  saturated  events,  (Bache,  et  al.,  1975),  whicn 
explains  their  coupling. 

3.4  HIGHER  ORDER  SOURCE  EFFECTS 

Many  potentially  important  source  effects  are  two-  or  three- 
dimensional  and  thus  cannot  be  represented  by  an  RDP  source.  These 
higher  order  effects  include  the  depth-dep*'-’"  ■-•vice  of  overburden 
pressure,  nonuniform  source  region  properties,  wlmear  interaction 
with  the  free  surface,  and  release  of  stored  strain  energy  near  the 
source.  Some  of  these  effects  can  be  approximated  with  elastic 
models  to  estimate  their  influence.  Results  from  these  kinds  of 
approximations  are  summarized  in  this  section. 

8ache,  Rodi  and  Mason  (1978)  used  elastic  models  to  attempt  to 
delineate  the  influence  of  higher  order  effects  including:  double¬ 
couple  generation  by  tectonic  release,  spallation,  attenuation  of 
upgoing  waves  by  spallation  and/or  scattering  in  the  near  surface 
materials.  This  leads  to  many  questions  about  the  physics  of  the 
source.  These  questions  can  probably  only  be  resolved  by  resorting 
to  multidimensional  calculations  like  those  we  describe  in  Section 
IV. 

Empirical  estimates  for  spall  closure  (e.g.,  Viecelli,  1973; 
Sobel,  1978)  suggest  that  the  generated  Rayleign  waves  might  be 
nearly  the  same  size  as  those  from  the  explosion  itself.  However, 
the  spall  induced  Rayleigh  wave  is  90  degrees  out  of  phase  with  that 
from  the  explosion  and  makes  an  unimportant  contribution  to  the  ALQ 
and  TUC  records  for  reasonable  assumptions  about  its  amplitude.  It 
is  not  clear  whether  this  result  can  be  generalized  to  allow  the 
spall  effect  to  be  ignored. 

A  potentially  important  spall  related  phenomenon  is  the  at¬ 
tenuation  of  upgoing  waves  from  the  source  during  the  spall 
process.  This  is  explored  by  Bache,  Rodi  and  Mason,  but  conclusions 
cannot  be  confidently  drawn  without  better  understanding  the  physics 
of  the  process.  Again,  two-dimensional  calculations  like  those 
described  in  Section  IV  are  needed  to  improve  our  understanding. 
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To  explore  the  effects  of  a  double-couple  component,  Bache, 
Rod!  and  Mason  (1978)  used  the  inferred  double-couple  solutions  of 
Toksoz  and  Kehrer  (1972).  In  this  form  the  double-couple  has  almost 
no  effect  on  the  waveform,  but  simply  scales  the  amplitude.  If  the 
events  studied  by  Toksoz  and  Kehrer  are  typical,  some  15  or  20 
percent  of  the  discrepancy  between  the  ALQ  and  TUC  inferred  RDP 
amplitudes  is  due  to  the  double-couple.  For  PILEORI VER,  their 
inferred  douole-couple  dominates  the  solution.  However,  their 
orientation  does  not  explain  the  ALQ-TUC  amplitude  discrepancy,  but 
actually  increases  it.  We  will  discuss  the  question  of  a  super¬ 
imposed  double-couple  in  more  detail  in  the  next  section. 

3.5  SURFACE  WAVE  RAOIATION  PATTERNS 
3.5.1  Introduction 

A  serious  problem  always  cited  to  discourage  too  much  reliance 
on  surface  waves  as  a  yield  indicator  is  the  observed  fact  that  many 
events  are  characterized  by  strong  surface  wave  radiation  patterns. 
The  conventional  wisdom  is  that  superimposed  on  the  explosion  source 
is  a  double-couple  aligned  with  the  local  stress  field.  From  a 
theoretical  point-of-view,  one  should  be  able  to  deduce  the  orien¬ 
tation  and  moment  of  this  double-couple  from  the  radiated  Love  and 
Rayleigh  waves  and  correct  the  data  for  this  effect.  This  was  the 
approach  followed  by  Bache,  Rodi  and  Mason  (1978),  described  in  the 
preceding  section.  They  used  the  double-couple  solutions  obtained 
by  Toksoz  and  Kehrer  (1972)  for  a  fairly  large  number  of  underground 
explosions. 

Toksoz  and  Kehrer  introduced  the  restrictive  assumptions  that 
the  double-couple  source  had  the  same  location  ano  time  function  as 
the  explosion  and  was  oriented  as  a  vertical  strike-slip  fault.  If 
enough  good  data  were  available,  these  assumptions  could  be 
relaxed.  For  example.  Rivers  and  von  Seggern  (1979)  have  inverted 
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for  the  best  moment  tensor  fit  to  the  PILEDRIVER  aata  using  the 

moment  tensor  inversion  scheme  of  Menaiguren  (1977).  Of  course,  any 

inversion  scheme  is  limited  by  the  accuracy  of  the  path  corrections 
used. 

How  well  does  the  explosion  plus  double-couple  model  fit  the 
data?  The  strike-slip  double-couple  usea  in  the  inversions  of 
Toksoz  and  Kehrer  does  not  give  a  particularly  good  fit.  Rivers  and 

von  Seggern  (1979)  seem  to  fit  the  PILEDRIVER  data  quite  well,  but 

with  a  model  for  which  the  physical  interpretation  is  not  entirely 
clear.  The  orientation  is  that  of  a  reverse  fault  and  both  the 
explosion  and  double-couple  are  at  a  depth  of  2.5  kilometers. 

It  does  seem  clear  that  the  explosion  plus  double-couple  model 
qualitatively  explains  the  long  period  data.  Perhaps  it  also  gives 
a  quantitative  explanation  that  can  be  used  to  correct  individual 
station  amplitudes.  8ut  more  detailed  work  is  necessary  to  prove 
this  concept.  We  have  looked  at  two  very  similar  events,  MIGHTY 
EPIC  ana  DIABLO  HAWK,  in  some  detail,  ana  the  results  illustrate 
that  the  data  are  not  easily  explained,  in  more  than  a  qualitative 
way,  by  superimposing  a  double-couple  on  the  explosion. 

3.5.2  Analysis  of  Regional  Motions  from  MIGHTY  EPIC  and  DIABLO  HAWK 

Bache,  Farrell  and  Lambert  (1979)  describe  an  analysis  of  the 
regional  ground  motions  of  the  Rainier  Mesa  explosions  DIABLO  HAWn 
ana  MIGHTY  EPIC.  These  were  very  similar  events  of  fairly  low 
yield.  The  recording  stations,  seven  for  MIGHTY  EPIC  and  ten  for 
DIABLO  HAWK,  were  broaa-band,  three-component  digitally  recoraea 
stations  operated  by  Lawrence  Livermore  Laboratory,  Sandia  Labor- 
atories  and  Systems,  Science  and  Software  (S  ).  The  distances 
were  131  to  368  kilometers. 

The  objective  of  the  experiment  was  to  infer  the  size  ana 
orientation  of  the  double-couple  component  associated  with  these 
events.  First,  the  data  were  processed  to  determine  the  following: 
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•  The  coherence  between  the  same  station  recordings 
of  the  two  events  is  very  high.  The  sources 
were,  therefore,  nearly  identical,  the  difference 
being  a  scale  factor,  which  may  be  azimuthal ly 
dependent. 

•  The  coherence  between  different  components  is 
low.  This  seems  to  indicate  that  the  horizontal 
component  waves  emanate  from  the  vicinity  of  the 
source. 

•  The  data  were  filtered  with  a  four  to  eight 
second  band  pass  filter  and  with  1*1SRS  to  identify 
the  fundamental  mode  Love  and  Rayleigh  waves  with 
a  period  near  six  seconds.  All  indications  are 
that  these  modes  were  successfully  isolated. 

The  best  fitting  explosion  RDP  plus  double-couple  was  deter¬ 
mined  by  fitting  the  LR,  IQ  and  LQ/LR  data.  The  double-couple  was 
assumed  to  be  a  vertical  strike-slip.  The  best  fitting  solutions 
for  DIABLO  HAWK  are  shown  in  Figures  3.7,  3.8  and  3.9.  The 
analogous  solutions  for  MIGHTY  EPIC  are  shown  in  Figures  3.10  and 
3.11. 

The  inferred  ROP,  (^J  is  within  the  range  of  expected  values 
for  such  events.  This  indicates  that  the  amplification  effects  of 
the  path  are  about  right  at  these  regional  distances.  The  fit  to  the 
Love  waves  is  not  very  good.  In  fact,  the  data  would  be  fit  nearly 
as  wel 1  with  a  circle! 

There  is  another  interesting  aspect  of  these  data.  The  DIABLO 
HAWK/MIGHTY  EPIC  ratio  is  about  0.9  at  high  frequencies  and  0.6  at 
low  frequencies.  This  is  true  for  nearly  all  components  at  all 
stations.  The  frequency  domain  ratios  can  also  be  seen  in  time 
domain  measurements  of  Pn  (ratio  of  0.9)  and  long  period  Love  and 
Rayleigh  waves  (ratio  of  0.6).  It  seems  that  the  total  (explosion 
plus  double-couple),  source  level  is  larger  at  long  periods  for 
DIABLO  HAWK.  There  is  no  easy  way  to  explain  this  in  the  context  of 
the  simple  theory  used. 
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Theoretical  Love  wave  radiation  patterns  are  compared  to  DIABLO  HAWK  observations 
normalised  to  2 50  km.  The  theoretical  solutions  have  6  =  90,  X  =  0,  strike  = 

NI6  E,  F  =  0.575,  ps  =  40  kbar  and  three  values  of  the  moment,  Mn. 


Figure  3.10.  The  fit  to  the  MIGHTY  EPIC  LQ/LR  data.  We  assume 
vertical  strike-slip  faulting. 
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3.5.3  Discussion  of  the  Radiation  Pattern  Data 

We  find  that  the  long  period  radiation  from  MIGHTY  EPIC  is 
about  twice  as  large  as  that  from  DIABLO  HAWK  on  all  three  com¬ 
ponents  at  nearly  every  station,  while  the  short  period  radiation  is 
about  the  same.  Also,  no  single  double-couple  source  can  provide  a 
very  good  fit  to  the  Love  wave  data  (which  are  almost  certainly  true 
Love  waves  originating  near  the  source,  though  the  amplitudes  may  be 
contaminated  by  off-azimuth  Rayleigh  wave  arrivals).  The  material 
properties  for  the  two  events  are  nearly  identical,  so  spectral  dif¬ 
ferences  in  the  explosion  source  functions  are  unlikely.  One  hypo¬ 
thesis  that  will  explain  these  data  is  that  the  MIGHTY  EPIC  event 
had  much  more  block  faulting  which  is  not  dominantly  along  one 
azimuth,  but  includes  slip  on  faults  at  several  azimuths.  If  the 
faulting  is  associated  with  tectonic  strain  release,  the  superim¬ 
posed  double-couple  sources  must  have  spectra  that  are  substantially 
smaller  at  1  Hz  than  at  longer  periods.  Since  the  double-couple 
source  dimensions  must  be  quite  small,  this  suggests  that  the  block 
faulting  must  occur  with  extraordinarily  slow  (compared  to  earth¬ 
quakes)  rupture  velocities,  on  the  order  of  100  m/sec. 

An  alternative  model  for  the  generation  of  shear  waves  that 
may  help  explain  the  observations  is  that  there  is  a  substantial 
amount  of  wave  conversion  near  the  source  due  to  passive  (i.e., 
releasing  little  or  no  stored  strain  energy)  fault  motions.  In  a 
recent  paper  Salvado  and  Minster  (1980)  have  shown  that  this  con¬ 
version  can  be  large  for  faults  that  are  weakly  bonded  and  close  to 
the  source.  However,  this  study  was  done  for  an  idealized  geometry, 
so  quantitative  estimations  of  the  potential  effect  on  the  seismic 
radiation  for  actual  events  has  not  yet  been  done. 

Part  of  the  difficulty  in  explaining  the  radiation  pattern 
with  a  model  for  tectonic  strain  release  is  that  the  models  used  are 
too  simple.  Stevens  (1980)  has  shown  that  an  unhomogeneous 
prestress  field  (including  stress  concentrations)  leads  to  a 
radiation  pattern  that  includes  higher  order  terms  than  the  double 
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couple.  One  conseouence  is  the  absence  of  nodes,  which  are  also 
absent  in  the  data.  This  is  another  more  realistic  effect  which  may 
help  explain  the  observed  complexity. 

3.6  CONCLUSIONS 

The  capability  exists  to  mot!  surface  waves  from  underground 
explosions  in  detail.  The  source  may  be  of  arbitrary  complexity  in 
plane  layered  earth  models.  In  this  section,  we  have  only  been  dis¬ 
cussing  elastic  point  source  models,  but  in  the  next  section  we  will 
show  that  complex  finite  difference  simulations  of  multidimensional 
explosion  sources  can  also  be  used  to  model  surface  waves.  An 
important  question  is  how  to  account  for  changes  in  the  local  source 
material  when  the  average  path  properties  clearly  remain  the  same. 
We  have  been  using  an  approximate  technique  that  seems  to  give 
reasonable  results. 

When  we  consider  the  source  amplitudes  inferred  from  ALQ  and 
TUC  data  using  these  detailed  modeling  procedures,  we  conclude  that: 

•  Even  with  tne  simplest  source  model  (RDP ) ,  the 
inferred  source  amplitudes  are  reasonably  con¬ 
sistent. 

•  For  similar  materials,  Rayleigh  wave  amplitude  is 

proportional  to  For  different  materials 

(e.g.,  tuff  and  granite),  the  dependence  on  u  is 
much  weaker  and  varies  with  period. 

•  Empirical  estimates  for  spall  closure  suggest 
that  the  generated  Rayleigh  waves  might  be  nearly 
the  same  size  as  those  from  the  explosion  it¬ 
self.  However,  spall  closure  appears  to  be 
unimportant  for  ALQ  and  TUC  records. 

Important  questions  that  require  further  investigation  include 
the  following: 


65 


SYSTEMS,  science  *na  SOFTWAIte 


Path  Effects 


•  What  are  the  effects  of  lateral  heterogeneities 
and  departure  from  plane  horizontal  layering? 
Are  such  effects  responsible  for  the  50  percent 
amplitude  difference  from  ALQ  to  TUC? 

•  Attenuation  (Q)  is  not  very  well-known,  though 
this  is  less  important  for  surface  waves  than 
body  waves. 

Source  Effects 

•  How  appropriate  is  the  one-dimensional  source 
representation?  Known  two- and  three-dimensional 
effects  include: 

1.  Depth  dependence  of  overburden  pres¬ 
sure. 

2.  Spallation  and  cracking  to  the  surface. 

3.  Block  motion,  induced  fault  motion, 
tectonic  strain  release. 

•  How  are  tangential  (SH  and  Love)  waves 

generated?  Why  are  they  generated  at  all 

azimuths?  What  parameters  control  their 
amplitude?  How  well  does  a  single  double-couple 
represent  this  perturbation  on  the  source? 

•  How  large  is  spall  closure?  How  are  the  free 

surface  reflected  waves  in  the  nonlinear 

environment  different  from  those  in  the  elastic 
approximation? 

•  The  ALO  and  TUC  study  raises  the  question,  if 

stations  at  many  azimuths  were  used,  how  much 
would  the  estimate  vary?  How  much  of  the 

"source"  radiation  pattern  is  a  path  effect?  How 
much  of  this  path  effect  can  be  corrected  by 
careful  modeling? 
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Finally,  we  conclude  that  a  great  deal  of  progress  has  been 
made  in  understanding  surface  wave  excitation  by  explosions.  There 
are  partial  answers  to  all  the  auestions  listed  above  and  the  tech¬ 
niques  are  available  to  obtain  more  complete  answers. 
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IV.  THE  CONTRIBUTION  OF  TWO-DIMENSIONAL  SOURCE  EFFECTS 
TO  THE  FAR-FIELO  SEISMIC  SIGNATURES  OF 
UNDERGROUND  NUCLEAR  EXPLOSIONS 


4.1  INTRODUCTION 

A  recent  report  by  Bache,  Barker,  Rimer  and  Cherry  (1980), 
which  has  the  same  title  as  this  section,  includes  an  analysis  of 
the  seismic  waves  from  a  series  of  axisymmetric  calculations  of 
explosions  in  granite.  The  abstract  of  that  report  is  reprinted 
below: 


Two-dimensional  calculations  of  underground  explo¬ 
sions  are  able  to  model  nonlinear  interaction  with  the 
free  surface  (including  spallation)  ana  the  depth  de¬ 
pendence  of  overburden  pressure  and  material  proper¬ 
ties,  as  well  as  the  complex  constitutive  behavior  of 
geologic  materials  which  has  been  extensively  studied 
with  one-dimensional  (spherically  symmetric)  calcu¬ 
lations.  In  this  study,  eleven  axisymmetric  calcula¬ 
tions  of  explosions  in  granite  are  analyzed  to  deter¬ 
mine  the  two-dimensional  effects  on  the  far-field  body 
and  surface  waves  ano  on  the  magnitudes  m^  and  Ms. 
Seven  of  these  are  hypothetical  150  KT  explosions  at 
depths  from  159  to  1000  meters  calculated  by  Applied 
Theory,  Inc.  (ATI).  These  show  depth  effects  that  are 
different,  though  not  strikingly  so,  from  those  cal¬ 
culated  with  one-dimensional  source  models.  The  other 
four  calculations  were  done  by  Systems,  Science  and 

Software  (S-Cu bed).  The  first  was  for  a  specific  event 
(PILEDRIVER) ,  ana  the  depth  and  yield  were  varied  for 
the  other  three.  The  near-  and  far-field  data  for 
PILEDRIVER  were  matched  rather  well,  except  that  the 
spallation  was  too  large  in  the  calculation. 

Synthetic  body  and  surface  waves  from  the  S-Cubed 
two-dimensional  sources  are  compared  to  those  from 

analogous  one-dimensional  sources  to  isolate  the  two- 
dimensional  effects.  These  are  not  important  for  the 
deep  events  or  for  the  first  arriving  P  wave  from  the 
shallow  events.  They  do  strongly  enhance  the  surface 
wave  amplitudes  (by  a  factor  of  two  or  three)  for  the 
two  shallow  calculations,  including  that  for  PILE- 

DRIVER,  and  reduce  the  m^  by  0.1  -  0.2  units.  The 

m*,  effects  can  clearly  be  attributed  to  spallation 
and  the  accompanying  reduction  of  the  pP  amplitude. 
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which  are  overpredicted  by  the  calculations.  The  Ms 
enhancement  must  be  due  to  the  same  causes,  though  they 
are  not  easily  isolated. 

The  first  section  of  the  Bache,  et  £l_.  (1980),  report  is  the 
"Introduction  and  Summary,"  which  induces  the  important  conclu¬ 
sions.  This  material  is  reproduced  in  the  following  subsections. 

4.2  BACKGROUND 

To  successfully  monitor  a  test  ban  treaty  prohibiting  or 
limiting  underground  nuclear  explosions,  it  is  necessary  to  unoer- 
stand  the  seismic  signatures  of  these  events.  An  important  part  of 
the  research  effort  to  improve  this  understanding  has  been  the 
development  and  application  of  deterministic  methods  to  compute  the 
seismic  wave  signatures  of  nuclear  explosions.  In  nearly  all  of 
this  work  the  source  is  assumed  to  be  spherically  symmetric  and  the 
earth  is  assumed  to  be  plane-layered.  Theoretical  seismograms  can 
then  be  computed  with  widely  available  methods  and  compared  to 
observations.  This  procedure  has  been  auite  successful  and  most  of 
the  important  controlling  parameters  have  been  identified  and  their 
effects  have  been  quantified  to  some  degree.  However,  many 
important  issues  remain  unresolved. 

Most  of  the  important  questions  that  remain  regarding  the 
generation  of  seismic  waves  by  underground  explosions  are  associated 
with  multi-dimensional  effects.  For  example,  there  is  not  yet  a 
clear  understanding  of  the  effect  on  the  seismic  radiation  of  nearpy 
interfaces,  the  free  surface  (allowing  spallation),  the  overburden 
pressure,  nonhydrostatic  prestress,  and  zones  of  weakness  in  the 
near  source  environment.  The  one-dimensional  calculations  now 
incorporate  detailed  constitutive  models  that  include  realistic 
models  for  pore  collapse,  effective  stress,  yielding,  and  cracking 
due  to  shear  and  tension  failure.  These  constitutive  models,  have 
been  generalized  to  two-dimensions  for  the  axisymmetric  finite 
difference  source  calculations  analyzed  in  this  report.  While  they 
do  not  include  all  the  multi-dimensional  effects  we  have  listed. 


they  do  include  some  of  tne  most  important,  ana  so  represent  a 
significant  step  forwara  in  tne  development  of  realistic  theoretical 
simulations  of  tne  seismic  waves  from  underground  explosions. 

4.3  AXISYMMETRIC  CALCULATIONS  OF  NUCLEAR  EXPLOSIONS 

In  this  report  we  present  a  detailed  analysis  of  the  seismic 
waves  from  eleven  two-dimensional  finite  aifference  calculations  of 
underground  nuclear  explosions  in  granite.  Seven  of  these  calcu¬ 
lations  were  done  by  J.  Trulio  and  N.  Perl  of  Applied  Theory,  Inc. 
(ATI)  and  four  by  N.  Rimer  and  J.  T.  Cherry  of  Systems,  Science  ana 
Software  (S-Cubed).  They  share  the  same  axisymmetric  geometry, 
thougn  the  details  are  different. 

Tne  two-aimensional ,  axisymmetric  explosion  calculations 
include  tne  presence  of  a  free  surface,  the  dependence  of  overburden 
pressure  on  depth  and,  for  the  S-Cuoed  calculations,  some  dependence 
of  material  properties  on  depth.  These  are  probaoly  tne  most 
important  higher  order  corrections  to  tne  one-dimensional  models. 
The  depth-dependence  of  overourden  pressure  is'  a  property  of  all 
test  sites.  Further,  our  understanding  of  geologic  structures  is 
normally  based  on  plane-layerea  models.  Therefore,  axisymmetry  is  a 
natural  geometry,  and  specification  of  tne  geometry  ana  material 
properties  in  two-dimensions  is  more  straightforward  than  char¬ 
acterizing  the  entire  near-source  environment  oy  a  heterogeneous 
material . 

Tne  most  important  two-dimensional  effects  are  associated  with 
the  nonlinear  interaction  of  the  stress  waves  witn  tne  free  sur¬ 
face.  Surface  spallation  is  an  obvious,  even  dominant,  phenomenon 
observed  in  the  near-field,  yet  it  has  never  been  included  in 

seismic  wave  propagation  studies  in  a  very  satisfactory  way. 

Perhaps  more  important,  tnere  is  ample  evidence  that  tne  f^ee 

surface  pnase  pP  is  more  complex  tnan  predicted  by  elastic  theory 

with  spherically  symmetric  sources,  but,  again,  this  remains  mostly 
in  tne  realm  of  speculation.  Witn  the  axisymmetric  calculations,  we 
are  aole  to  study  these  important  effects. 
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The  objectives  of  the  parameter  variations  in  the  ATI  and 
S-Cubed  suites  of  calculations  are  somewhat  different.  The  ATI 
calculations  were  all  done  at  the  same  yield,  150  kt,  in  a 
hypothetical  granite  halfspace.  Only  the  depth  was  varied,  from  159 
to  1000  meters.  The  philosophy  of  the  S-Cubed  calculations  was  to 
begin  by  modeling  a  specific  event,  P ILEDR I VER ,  at  the  Nevada  Test 
Site.  The  constitutive  model  and  source  geology  (three  layers)  were 
chosen  for  this  event.  The  computed  and  observed  ground  motions 
were  compared  at  some  twenty-five  near-field  gauge  locations.  The 
comparison  was  quite  good,  except  that  the  calculation  overpredicteo 
the  amount  of  cracking  (spallation)  near  the  surface.  This  turns 
out  to  be  important  when  comparing  synthetic  and  observed  short 
period  data. 

The  other  three  S-Cubed  source  calculations  were  the  same, 
except  the  depth  and  yield  were  varied  from  the  values  appropriate 
for  PILEORIVER.  The  four  calculations  were  as  follows: 

Depth  (Meters)  Yield  (Kt) 


463 

1000 

1000 

400 


60  (PILEDRIVER) 
150 
20 
20 


4.4  OUTLINE  OF  THE  ANALYSIS 

The  analyses  of  the  seven  ATI  and  four  S-Cubed  calculations 
are  described  separately,  but  follow  parallel  lines.  First,  we 
describe  the  near-field  ground  motions  predicted  by  the  calculations 
in  some  detail.  Then  we  present  theoretical  seismograms  for  these 
calculations  and  analyze  their  implications.  This  outline  is 
apparent  in  the  section  headings  listed  below: 
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Section  II: 


ATI  Granite  Calculations 


Section  III:  Ms  ana  Estimates  for  ATI 

Granite  Calculations 

Section  IV:  S-Cubea  Granite  Calculations 

Section  V:  Ms  and  015  Estimates  for  S-Cubed 

Granite  Calculations 

Section  VI:  Comparison  of  Predicted  and  Observed 

Seismograms  for  P ILEDR I VER 

Section  VII:  Comparative  Analysis  of  ATI  and 
S-Cubed  Granite  Calculations. 

A  key  step  in  this  study  is  the  linkage  between  the  near¬ 
field  ground  motions  computed  by  the  finite  difference  programs  ana 
the  analytical  techniques  used  to  propagate  seismic  waves  in 
realistic  earth  models.  A  rigorous  method  for  accomplishing  this 
linkage  is  described  in  Appendix  A,  "Synthetic  Seismograms  from 
Complex  Source  Calculations."  The  details  and  examples  presentea  in 
this  appendix  are  mostly  for  the  calculation  of  the  normal  mode 
(Rayleigh  waves)  displacements  in  a  plane-layered  earth  model,  but 
ray  theory  methods  for  propagating  body  waves  can  be  used  within  the 
same  theoretical  framework. 

To  compute  synthetic  seismograms  at  large  distances  for  the 
finite  difference  source  calculations,  it  is  necessary  to  monitor 
the  tractions  and  oisplacements  on  some  (hypothetical)  surface  which 
entirely  encloses  the  region  of  nonlinear  material  response.  Our 
discussion  of  the  source  calculations  in  Sections  II  ana  IV  is 
concerned  with  the  characteristics  of  the  monitored  ground  motions 
on  this  elastic  surface. 

An  important  constraint  on  the  numerical  results  is  the  re¬ 
quirement,  based  on  conservation  of  momentum,  that  the  total  down¬ 
ward  force  and  impulse  (on  the  surface  enclosing  the  source  region) 
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vanish  at  late  time.  It  turns  out  (Appendix  A)  that  the  inevitable 
numerical  errors  that  cause  deviation  from  this  requirement  can 
dominate  the  solution  at  long  periods.  Therefore,  a  major  theme  in 
Sections  II  and  IV  is  the  application  of  a  "correction"  to  adjust 
the  computed  vertical  force  and  impulse  to  zero  at  the  last  time 
step. 

The  synthetic  seismogram  results  are  presented  in  Sections  III 
and  V.  The  vertical  force  "correction"  is  also  an  issue  in  these 
sections,  because  it  would  be  unsatisfactory  for  an  ad_  hoc  cor¬ 
rection,  which  is  what  we  apply,  to  dominate  the  answer.  Our 
conclusion  is  that  the  vertical  force  terms,  after  correction,  do 
not  play  an  important  role.  All  of  the  solution  checks  we  are  able 
to  make  indicate  that  the  body  and  surface  wave  synthetic 
seismograms  accurately  represent  the  seismic  waves  generated  by 
these  theoretical  sources. 

The  synthetic  seismogram  analyses  in  Sections  III  and  V  are 
focussed  on  comparing  the  seismic  waves  from  the  two-dimensional 
sources  with  those  from  analogous  one-oimensional ,  reduced  displace¬ 
ment  potential  (RDP)  sources.  For  the  S-Cubed  calculations,  we  are 
able  to  compare  to  the  RDP  source  computed  with  the  same  constitu¬ 
tive  model  in  spherical  symmetry.  No  such  one-dimensional  calcu¬ 
lation  was  available  for  the  ATI  granite,  so  we  compare  to  the  RDP 
predicted  by  the  semi-empirical  model  of  Mueller  and  Murphy  (1971). 

Most  of  the  calculations  are  for  hypothetical  events,  so 
direct  comparison  with  observed  seismograms  is  only  possible  for  the 
PILEDRIVER  calculation  done  by  S-Cubed.  This  comparison  is  made  in 
Section  VI.  Finally,  in  Section  VII,  the  results  of  all  eleven 
calculations  are  plotted  together  for  direct  comparison. 

4.5  CONCLUSIONS 

This  study  is  primarily  an  investigation  of  the  influence  of 
burial  depth  on  the  seismic  signals  from  underground  explosions. 
For  the  most  part,  the  results  are  interesting,  but  not  terribly 


exciting,  because  they  are  pretty  much  in  accord  with  our  expec¬ 
tations.  The  key  exception  is  the  surface  wave  results  for  tne 
S-Cubed  calculations. 

The  most  important  result  of  this  study  is  that  the  S-Cubea 
calculations  show  M$  to  be  a  strong  function  of  depth.  The 
shallow  S-Cubed  calculations  have  surface  wave  amplitudes  that  are  a 
factor  of  two  or  three  larger  than  those  from  a  comparable  one-di¬ 
mensional  source  calculation.  This  effect  is  probably  exaggerated 
because  the  free  surface  interaction  effects  are  too  large  in  the 
S-Cubed  calculations,  but  would  remain  important  even  if  the  near 
surface  material  were  strengthened.  On  the  other  hand,  the  ATI 
calculations  show  no  strong  dependence  of  on  depth,  even  though 
the  shallow  ATI  sources  cratered.  We  do  not  know  why  the  two  sets 
of  calculations  give  such  different  results. 

The  important  conclusions  are  listed  at  the  end  of  several 
sections  of  the  report.  Summaries  for  the  ATI  and  S-Cubea  calcu¬ 
lation?  are  given  separately  in  Sections  3.7  and  5.10.  Our  conclu¬ 
sions  about  the  comparison  of  synthetic  and  observed  seismograms  are 
listed  in  Section  6.4,  and  all  of  Section  VII  shoula  be  read  as  a 
summary. 

The  main  results  are  listed  below: 

ATI  Calculations 

•  Neither  nor  are  strongly  aependent  on 
depth.  The  most  effect  was  on  M$  at  shallow 
depths. 

Compared  to  the  Muel  ler /Murphy  RDP  source,  the 
depth  dependence  in  both  amplitude  and  corner 
frequency  is  less  for  the  ATI  sources. 

•  The  pP  phase  appears  to  be  smaller  than 
expected  from  elastic  theory. 
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S-Cubed  Calculations 


•  Two-dimensional  effects  are  not  very  important 
for  the  two  deep  explosions  (20  kt  and  150  kt 
at  1000  meters).  Both  m^  and  M$  are  little 
different  from  the  values  estimated  from  an  ROP 
source  computed  with  the  same  constitutive 
model  at  the  same  depth.*  This  is  true  even 
though  considerable  cracking  and  spallation 
occur  in  the  two-dimensional  calculations. 

•  For  body  waves  the  first  arriving  P  wave  is 
essentially  the  same  for  one-  and  two-di¬ 
mensional  sources.* 

•  The  two-dimensional  effects  enhance  the  surface 
wave  amplitudes  for  the  shallow  events  (60  kt 
at  463  meters  and  20  kt  at  400  meters)  by  a 
factor  of  two  or  three.  We  must  qualify  this 
by  pointing  out  that  these  shallow  calculations 
have  very  strong  surface  interaction  effects. 

Comparison  with  PILEORIVER  data  indicates  that 
the  free  surface  interaction  is  overpredicted, 
at  least  for  that  event.  Less  free  surface 
interaction  would  presumably  give  less 
enhancement  of  the  surface  wave  amplitudes. 

•  For  the  shallow  events,  the  mb  is  different 
than  predicted  with  an  RDP  source,  though  by 
less  than  0.2  units. 

•  Analysis  of  the  spectra  show  that  in  no  case  is 
pP  a  spectral  shadow  of  P,  as  it  is  for  an  ROP 
source  and  elastic  propagation. 

*  these  results  are  essentially  tests  of  the  entire  computational 
procedure.  It  is  gratifying  that  these  two  very  different  and 
complex  procedures  arrive  at  the  same  results. 
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•  Phases  that  seem  to  be  associated  with  spall 
closure  can  be  seen  on  the  shallow  source  body 
wave  records.  However,  they  are  not  easily 
associated  with  identifiable  crack  closure 
patterns  in  the  source  calculation. 

The  comparison  of  observed  and  calculated  body  and  surface 
waves  for  PILEDRIVER  in  Section  VI  leads  to  the  conclusion  that  the 
two  are  in  rather  good  agreement.  This  conclusion  is  subject  to  the 
Qualifications  one  often  faces  in  this  kind  of  comparison.  For  the 
surface  waves,  it  is  the  need  to  account  for  the  non-ax i symmetric 
component,  usually  attributed  to  tectonic  stress  release.  Adding  a 
recent  estimate  for  this  component  by  Rivers  and  von  Seggern  (1979) 
to  our  solution,  we  get  good  agreement  with  the  data. 

For  body  waves  the  comparison  is  complicated  by  the  apparent 
presence  of  strong  azimuthal  effects  in  the  radiated  short  period 
energy  (Hadley  and  Hart,  1979).  However,  our  conclusion  is  that  the 
computed  PILEDRIVER  source  (in  one-  or  two-dimensions)  has  about  the 
right  direct  P  amplitude.  However,  the  two-dimensional  source 
calculation  appears  to  include  too  much  non-linear  interaction  with 
the  free  surface.  The  lat,er  portion  of  the  P  waveform  does  not 
match  the  data,  apparently  because  pP  is  too  greatly  suppressed  ano 
because  the  seismic  energy  from  spall  closure  is  too  large  or  is 
timed  incorrectly.  This  "overprediction"  of  surface  interaction 
effects  is  expected  since  comparison  of  theoretical  and  ooserved 
near-field  motions  and  plots  of  the  cracking  near  the  source 
indicate  that  there  was  too  much  spallation  in  the  calculation. 

All  the  theoretical  and  mb  values  are  plotted  together 
in  Section  VII.  They  are  shown  versus  source  depth  and  versus 
yield,  including  the  observed  values  for  HARDHAT,  SHOAL  and  PILE- 
DRIVER.  The  most  dramatic  difference  between  the  ATI  and  S-Cubeo 
calculations  is  the  strong  dependence  of  M$  on  depth  predicted  by 
S-Cubed,  but  not  by  ATI  sources.  This  must  be  a  reflection  of  the 
different  constitutive  models  used. 
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The  constitutive  models  used  by  S-Cubed  (in  spherically 
symmetric  source  calculations)  lead  to  RDP  source  functions  that  are 
strongly  peaked,  with  the  value  near  1  Hz  a  factor  of  five  or  more 
larger  than  the  value  at  long  periods.  The  peaking  is  due  to  the 
incorporation  of  an  effective  stress  law  and  the  choice  of 
unconfined  compressive  strength  (0.75  kbar),  based  on  laboratory 
data  for  fractured  granite  and  results  of  comparison  with  near-field 
ground  motion  observations  (Rimer,  personal  communication). 

While  we  do  not  have  RDP  source  functions  for  the  ATI  granite, 
comparison  of  M$  and  m^  for  the  ATI  two-aimensional  calculations 
indicates  that  the  RDP  peaking  is  probably  less  than  a  factor  of 
two.  Due  primarily  to  this  difference,  the  m^  for  the  S-Cubed 
calculations  is  about  0.5  units  higher  than  that  for  ATI  calcula¬ 
tions  of  the  same  yield.  The  ATI  values  fall  between  those  for 
the  shallow  and  deep  S-Cubed  calculations. 


V.  THREE-DIMENSIONAL  EARTHQUAKE  MODELING  INCLUDING 
NONLINEAR  RUPTURE  DYNAMICS 


We  are  using  a  th ree-a i men si on a 1  finite  difference  method  to 
model  those  aspects  of  earthquake  dynamics  believed  to  be  important 
for  predicting  the  teleseismic  signal.  The  calculations  nave  simu¬ 
lated  inelastic  response  of  the  fault  zone  through  nonlinear  bound¬ 
ary  conditions  on  the  fault  plane.  The  numerical  method  admits  more 
general  nonlinear  material  behavior,  which  is  being  incorporated  as 
required. 

Our  initial  numerical  simulations  were  for  a  simple,  constant 
rupture  velocity,  constant  stress-drop  earthquake  mooel.  These 
results  show  a  significant  influence  of  fault  width  on  the  slip 
function.  For  points  on  the  fault  which  are  more  than  a  fault  width 
away  from  the  focus,  the  slip  function  is  spatially  quite  uniform. 
The  static  slip  is  given  by  aow/pa1-,  in  agreement  with  estimates 
based  on  elastostatic  solutions.  The  numerical  simulations  predict 
a  rise  time  of  w/2a  and  peak  slip  velocity  of  >/2fw/B'  ao/ps.  In 
these  expressions,  w  is  the  fault  width,  a<r  the  stress-drop,  f  the 
maximum  frequency,  s  the  shear  speed,  and  o  the  density. 

Careful  studies  of  the  broad-band  character  of  earthquake 
teleseismic  signatures  indicate  that  the  relative  excitation  of 
long-  and  short-period  waves  is  inconsistent  witn  predictions  from 
the  simple  model.  In  addition,  a  few  isolated,  severe  stopping 
phases  dominate  near-field  accelerograms  synthesized  from  the  simple 
model,  in  conflict  witn  strong  motion  data.  More  complex  rupture 
physics  is  required  to  explain  the  discrepancies  between  the  model 
predictions  and  teleseismic  and  strong  motion  observations.  In 
particular,  a  better  physical  description  of  the  cessation  of 
rupture  growth  appears  necessary. 

A  displacement-weakening  constitutive  model,  described  in 
Figure  5.1,  has  been  applied  to  simulate  spontaneous  rupture  in  the 
presence  of  both  uniform  and  heterogeneous  prestress.  In  this  model 
of  failure,  the  shear  strength  of  the  fault  zone  is  assumed  to  be  a 
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decreasing  function  of  its  deformation;  energy  is  dissipated  during 
rupture  extension  and  shear  stress  is  hounded  on  the  fault  edges. 
When  rupture  initiates  in  a  limited  region  of  high  stress  (Figure 
5.2a),  the  model  leads  to  spontaneous  propagation,  deceleration,  and 
stopping  of  faulting.  Figure  5.3  shows  the  evolution  of  the  rupture 
front  in  this  case,  as  well  as  the  far-field  P  and  S  wave  oisplace- 
ment  spectra  and  the  corresponding  displacement,  velocity,  and 
acceleration  pulses.  In  the  presence  of  multiple  hign-stress 
regions  (Figure  5.2b),  each  isolated  stress  concentration  behaves 
roughly  as  an  independent  event.  As  shown  in  Figure  5.4,  rupture 
velocity  becomes  quite  irregular,  and  is  strongly  coupled  to  peak 
velocity  of  slip  on  the  fault  surface.  Figure  5.5  shows  the 
irregular  rupture  growth,  as  well  as  the  far-field  spectra  and 
time-domain  pulses  which  result.  Predicted  accelerograms  are  in 
better  qualitative  accord  with  strong  motion  data. 

The  model  predicts  significant  consequences  for  the  spectral 
and  time-domain  characteristics  of  the  earthquake  teleseismic  signal 
due  to  the  presence  of  stress  irregularities.  For  example.  Figure 
5.6  compares  the  short-period  teleseismic  P-wave  for  the 
stress-concentration  model  with  that  for  a  uniform  stress  model 
(scaled  to  a  source  dimension  equal  to  the  dimension  of  the  stress 
concentration) .  Both  the  period  and  amplitude  of  the  b  phase  (first 
peak-to-trough)  agree  for  the  two  models,  while  the  stress-con¬ 
centration  model  has  twice  the  moment  of  the  uniform  stress  model. 
The  two  events  have  been  superimposed  on  an  mb  versus  plot 
(from  Filson  and  Bungum,  1972)  to  illustrate  that  the  stress-concen¬ 
tration  phenomenon  enhances  separation  of  earthquake  and  explosion 
populations. 
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Figure  5.R  Stress-concentration  fault  models  considered  in  this 
study. 
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Figure  S.6  Short-period  P-waves  from  uniform  stress  model  and  stress-con¬ 
centration  model.  The  two  events  are  plotted  as  triangles  on  an 
m|}-versus  Ms  diagram  (from  tilson  and  Bungum),  in  which  open 
circles  and  closed  circles  represent  presumed  explosions  and 
earthquakes,  respectively. 
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VI.  AUTOMATED  MAGNITUDES,  mb  AND  M$ 

6.1  INTRODUCTION 

In  this  section  we  describe  algorithms  for  computing  automated 

A  A 

spectral  body  and  surface  wave  magnitudes,  mb  and  M$,  that 
promise  to  give  more  consistent  and  convenient  measures  of  signal 
energy  than  the  currently  used  time  domain  magnitudes.  The  use  of 
these  magnitudes  is  illustrated  by  application  to  seismograms 
constructed  by  embedding  synthetic  signals  in  actual  noise 
recordings. 

The  automatic  magnitudes  may  be  viewed  as  a  byproduct  of 

3 

signal  processing  by  the  S  MARS  program.  This  program  is  based 
on  the  application  of  a  series  of  Gaussian  narrow-bana  filters  to 
the  data.  Applications  include  the  following: 

1.  Determination  of  phase  and  group  velocity 
dispersion  of  surface  waves.  An  example  of  the 
application  of  MARS  for  this  purpose  was 
mentioned  in  Section  3.3. 

2.  Detection.  MARS  was  implemented  as  a  P  wave 
detector  during  the  VSC  conducted  discrim¬ 
ination  experiment  in  FY  '78  and  '79. 

3.  Discrimination.  The  MARS  program  computes  high 
and  low  frequency  spectral  estimates  called 

3 

mb(f).  The  discriminant  used  by  S  in  the 
discrimination  experiment  is  based  on  compari¬ 
son  of  these  mb(f)  values  with  earthquakes 
and  explosions  falling  in  different  portions  of 
the  plane. 

A  natural  extension  of  this  work  is  to  use  MARS  to  auto¬ 
matically  provide  the  magnitude  needed  to  estimate  the  size  of  the 
event.  Ultimately,  the  program  could  automatically  detect,  discrim¬ 
inate  and  estimate  yield. 


i 
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Bache  (1979)  and  Bache,  Oay  and  Savino  (1979)  proposed  an 
algorithm  for  determining  a  MARS  based  body  wave  magnitude  called 
mb.  This  algorithm  was  tested  by  processing  recordings  of  eleven 
Pahute  Mesa  explosions  from  six  teleseismic  stations.  The  mb 
values  were  compared  to  a  carefully  determined  time  aomain  m^. 

The  statistics  of  a  linear  regression  on  log  yielo  were  also 
compared.  We  concluded  that  the  m^  was  at  least  as  gooa  a 
magnitude  measure  and  yield  indicator  as  the  most  carefully 
determined  time  domain  mb  for  the  high  signal/noise  data  processed. 

The  m^  algorithm  used  in  the  earlier  reports  worked  very 
well,  but  was  primitive  because  it  failed  to  account  for  the 

presence  of  seismic  noise.  The  discrimination  experiment  work  led 
to  the  development  of  more  sophisticated  algorithms  for  using  MARS 
output  to  determine  spectral  amplitude  (Masso,  et  a]_.,  1979;  Savino, 
et  £l_. ,  1980).  These  algorithms  incorporate  corrections  for  the 
presence  of  interfering  phases  and  seismic  noise. 

In  Section  6.2  we  describe  a  more  advanced  algorithm  for  com¬ 
puting  mb<  The  properties  of  this  new  mb  should  be  nearly  the 
same  as  the  properties  of  the  mb  used  in  the  earlier  studies.  The 
algorithm  is  tasted  by  applying  it  to  compute  mb  for  a  synthetic 
seismogram  embedded  in  varying  levels  of  seismic  noise.  Accurate 

results  are  obtained,  even  when  the  S/N  ratio  is  so  low  that  a  time 
domain  mb  measurement  would  probably  not  be  attempted.  We  com¬ 

plete  our  discussion  of  mb  by  computing  it  for  a  number  of  obser¬ 
vations  of  Soviet  explosions  recorded  at  RKON. 

A 

In  Section  6.5  we  describe  the  M$  algorithm,  which  is  much 
the  same  as  the  mb  algorithm,  except  that  the  dispersion  character 
of  the  surface  wave  signal  of  interest  must  be  specified  explic¬ 
itly.  The  time  series  is  first  processed  so  a  signal  with  the 

A 

expected  dispersion  appears  to  be  nearly  undisperseo,  and  the  M$ 
is  then  computed  just  like  mb< 

A 

Tests  of  the  M  algorithm  are  discussed  in  Section  6.6.  As 

A  *  A 

with  mb>  these  involve  computation  of  M$  for  a  synthetic  seismo¬ 
gram  embedded  in  varying  levels  of  seismic  noise.  We  summarize  our 
conclusions  in  Section  6.7. 
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6.2  THE  mb  ALGORITHM 

The  basic  algorithm  for  computing  mb  is  the  same  as  that 
used  by  Savino,  et  al_.,  (1980)  to  compute  the  high  and  low  frequency 
m^f)  used  for  discrimination.  In  fact,  it  is  nearly  correct  to 
say  that  mb  *  mb  (1  Hz).  There  are  some  details  that  are  dif¬ 
ferent,  as  we  shall  point  out  when  outlining  the  algorithm.  The 
steps  are  as  follows: 

1.  The  seismogram  is  processed  by  MARS,  which 
applies  a  suite  of  narrow-band  filters  to  the 
Fourier  transformed  seismograms.  The  peaks  of 
the  envelope  functions  are  saved  for  further 
processing. 

2.  The  MARS  detection  algorithm  is  applied  to 
identify  an  undispersed  P  wave  arrival  The 
exact  form  of  the  detection  algorithm  is  a 
subject  for  current  research,  but  we  used  the 
particular  algorithm  chosen  by  Savino,  et  aK , 

(1980),  which  will  be  discussed  in  more  detail 
below. 

3.  A  particular  frequency  is  chosen  for  computing 
a  "variable  frequency  magnitude,"  or  mb(f). 

We  chose  1  Hz,  though  this  may  not  be  the  best 
choice.  We  will  describe  the  details  of  the 
mb  (1  Hz)  calculation  in  more  detail  below. 


Detection 

The  detection  algorithm  used  by  Savino,  et  a1_. ,  (1980)  is 
illustrated  in  Figure  6.1.  The  envelope  function  peaks  are  plotted 
at  the  associated  group  arrival  time  (t^)  for  each  filter  center 
frequency  (fc).  The  amplitudes  (Ac)  are  indicated  by  the  symbol 
on  the  plot,  with  an  asterisk  denoting  the  largest  peak  for  each 
center  frequency.  The  other  symbols  are  0  to  9  with,  for  example,  a 
7  indicating  a  peak  with  amplitude  between  70  and  80  percent  of  the 
largest  peak  for  that  center  frequency. 
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Figure  6.1.  Bottom  portion  Is  a  seismogram  recorded  at  a  station 
In  Bluff,  Alaska,  from  an  event  in  the  Kurils.  The 
top  half  is  a  plot  of  filter  center  frequency  versus 
group  arrival  time  (sec).  The  t*  with  arrow  attached 
shows  where  the  automated  detection  and  first  arrival 
algorithms  picked  the  signal  arrival  time. 
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The  detection  algorithm  starts  with  a  "surface1*  of  the  kind 

shown  in  the  figure  and  searches  for  a  "ridge"  which  indicates  an 

undispersed  arrival.  This  pattern  recognition  is  done  iteratively, 

with  the  time  window  being  progressi vely  narrowed  about  regions  that 

appear  to  include  a  signal  witn  the  desired  characteristics.  The 

t  -f  section  shown  in  Figure  6.1  is  the  final  pattern  used  for 
9  c 

defining  the  group  arrival  time,  t^,  of  the  detected  signal.  Tnis 
is  the  time  of  arrival  of  the  maximum  signal  energy.  The  algorithm 
also  identifies  the  first  arrival  time  of  the  signal,  t*  ,  which 
is  noted  on  the  seismogram. 

Spectral  Resolution 

To  compute  the  variable  frequency  magnitude  values,  it  is 
necessary  to  identify  one  spectral  amplitude  (Ac)  for  each  f  . 
In  figure  6.1  tne  and  signs  denote  an  "acceptance  window" 
centered  about  the  group  arrival  time  t^.  The  window  is  based  on 
some  statistical  measures  of  the  uncertainty  in  the  t  estimate. 
Any  Ac  within  this  window  could  be  associated  with  the  detected 
signal,  and  there  are  often  several  amplitudes  for  a  particular  f 

V 

Savino,  et  aK ,  (1980)  select  one  Ac  for  each  fc  in  the 

following  way.  Each  amplitude  within  the  acceptance  window  is  given 
a  score  based  on  its  amplitude  and  time  separation  from  tg.  The 
weighting  factors  were  chosen  so  the  largest  Ac  usually  wins.  The 
time  separation  plays  a  role  in  selecting  between  points  of  similar 
ampl itude. 

This  spectral  resolution  algorithm  is  somewhat  arbitrary  and  a 
better  algorithm  might  occasionally  choose  a  different  Ac  for  a 
particular  signal.  One  weakness  that  becomes  apparent  in  the  mQ 
calculations  we  will  be  discussing  is  that  tne  algorithm  chooses  an 
Ac  for  each  f£  independent  of  the  choice  at  nearby  f  .  The 
algorithm  used  by  Bache  (1979)  also  used  a  score  to  choose  the  A  , 
but  the  score  was  intended  to  select  points  tnat  lead  to  the 
"smoothest"  dependence  of  Ac  and  t^  on  fc.  Of  course,  all 
reasonable  algorithms  give  the  same  choice  of  A  (f  )  when 
signal /noise  is  large. 
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Variable  Frequency  Magnitude  Calculation 

Having  identified  one  A  for  each  f  ,  we  have  an  estimate 

t  v* 

for  the  spectrum  of  the  detected  signal.  For  discrimination, 
Savino,  et  aU  (1980),  compute  m^(f)  from  these  Ac  at  two  fre¬ 
quencies,  one  high  (e.g.,  2.0  Hz)  and  one  low  (e.g.,  0.5  Hz).  The 
magnitude  is  computed  from 

mb(fc)  «  log  (A**(fc)  *  fc)  +  8,  (1) 

where  8  is  the  usual  distance  correction.  The  A^**  is  the 
envelope  amplitude,  Ac>  corrected  for  noise  (which  can  also  be 
nearby  "signals").  This  is  done  by  subtracting  the  contribution  of 
nearby  envelope  peaks  from  the  Ac  for  the  selected  peak.  Also, 
rather  than  use  the  noise-corrected  amplitudes  at  a  single 
freauency,  a  band  of  frequencies  is  selected.  A  noise-weighted, 
least  squares  quartic  fit  to  the  amplitudes  in  that  band  is  then 
made,  ano  the  Ag**(fc)  is  the  value  of  the  fittea  quartic  at 
fc.  Thus,  a  kind  of  “smoothing"  is  applied  to  tne  spectrum. 

Noise  Uncertainty 

We  mentioned  the  "deterministic"  noise  correction  that  ac¬ 
counts  for  the  influence  of  nearby  envelope  peaks  on  the  amplitude 
of  the  selected  envelope  peak.  This  correction  accounts  for  nearby 
"signal"  as  well  as  true  seismic  noise.  Random  noise  is  also 
important,  so  the  noise  spectrum  is  estimated  from  a  noise  winaow, 
usually  50  seconds  or  so,  preceding  the  detected  signal.  The  larger 
the  ambient  noise,  the  larger  the  uncertainty  in  the  amplitude 

A  **(f„).  A  statistical  estimate  of  this  uncertainty  is  made  as 
g  c 

part  of  the  m^(f)  calculations. 
mK  Calculation 

To  compute  the  m^  of  this  report,  we  essentially  calculated 
3^(1  Hz).  The  calculation  is  illustrated  in  Figure  6.2.  The 
seismogram  is  a  synthetic  for  an  explosion  source.  It  is  added  to  a 
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FILTER  CENTER  FREQUENCY  (Hz) 


Figure  6.2.  Calculation  of  for  a  seismogram  with  peak  synthetic/ 
peak  noise  *  100.  The  seismogram  is  shown  at  the  top 
with  an  asterisk  denoting  the  signal  arrival  time  t_. 


noise  sample  taken  from  an  RKQN  recording  in  tne  AI  data  base 
(Savino,  et  al^.,  1980).  For  this  case  the  ratio  of  the  largest  peak 
on  the  synthetic  to  the  largest  peak  on  the  noise  sample  (peak  syn- 
thetic/peak  noise)  is  100. 

The  mu  calculation  is  shown  graphically  at  the  bottom.  The 

D 

heaviest  line  is  a  plot  of  log  (Ag**  *  fc)  +  B  versus  f_.  Tne 

t  associated  with  each  A  **  is  plotted  as  the  lightest  line, 
9  9  = 

with  the  scale  at  the  right.  Zero  on  the  tg  scale  is  tg,  which 

is  also  indicated  by  an  asterisk  on  the  seismogram.  The  lowest  line 

of  medium  shade  is  log( AN  •  fc)  ♦  B,  where  A^  is  the  peaic 

noise  amplitude  at  fc» 

The  main  information  about  is  printed  with  the  graph. 

The  mu  is  cal  lea  MB  and  is  marked  with  an  X  on  the  spectral  plot. 
0 

The  noise  uncertainty  is  *  DMB.  These  values  are  not  symmetric 

about  m^,  though  they  appear  to  be  when  they  are  small  as  in  this 
example.  The  B  correction  used  to  convert  to  magnituae  is  listed  as 
BDEl. 

The  m.  in  this  and  subseauent  examples  in  this  report  is 
0 

computed  from  the  noise  corrected  amplitude  at  1  Hertz.  Men 
describing  the  mb(f)  algorithm,  we  mentioned  that  the  amplitude  at 
the  selected  frequency  is  taken  from  a  least  sauares  ouartic  fit  to 
the  amplitudes  in  a  band  about  that  freauency.  Since  we  used  only 
five  points  in  the  band,  this  part  of  the  algorithm  is  inoperative 
for  our  mb  examples.  The  intent  was  to  smooth  the  spectrum  over 
the  band  from  0.8  to  1.2  Hertz,  marked  by  vertical  bars  in  the  graph 
and  listed  as  FLEFT  and  FRIGHT.  In  the  examples  to  be  shown  a 
smooth  fit  can  be  made  by  eye  and  we  can  see  that  it  would  make 
little  difference  in  almost  every  case. 

The  remaining  undefined  quantities  with  the  mb  graph  are 
MBNF  and  S/N.  The  M8NF  is  the  noise  magnitude  at  1  Hertz.  The  S/N 
is  the  mean  signal/noise  power  over  the  selected  band  (0.8  to  1.2 
Hz). 

The  mb  for  the  seismogram  in  Figure  6.2  is  3.589.  This  may 
be  compared  to  the  time  domain  mb  of  3.776,  which  is  computed  from 
log  (A/T)  +  8,  with  T  -  0.72  seconds.  Systematic  differences 

93 


SYSTEMS.  SCIENCE  AND  SOFTWARE 


between  the  time  domain  m^  and  the  spectral  m^  are  expected.  We 
believe  the  m^  is  a  more  convenient  and  consistent  measure  of  the 
spectral  energy  in  the  P  wave. 

The  Ag**  spectrum  from  which  the  m^  is  determined  is 
simply  a  smoothed  version  of  the  Fourier  spectrum  of  an  isolated 
pulse.  This  is  illustrated  in  Figure  6.3,  where  the  Fourier  spec¬ 
trum  of  the  original  synthetic  seismogram,  which  was  embedded  in 
noise  to  make  the  seismogram  in  Figure  6.2,  is  plotted  with  the  MARS 
constructed  spectrum  from  that  figure. 

We  have  described  the  m^  algorithm,  as  currently  imple¬ 
mented.  There  are  several  aspects  that  could  be  improved.  In 
particular,  the  amplitude  used  in  the  magnitude  calculations  should 
be  based  on  the  fit  to  values  in  a  band,  rather  than  a  single  fre¬ 
quency  value.  Also,  the  algorithm  for  selecting  the  appropriate 
amplitudes  for  each  f  might  be  improved  by  including  some 
information  about  the  "smoothness"  of  the  frequency  dependence. 

In  the  next  section  we  demonstrate  that  this  algorithm  gives 
excellent  results  in  its  current  form.  Thus,  we  cannot  expect  large 
improvements  from  any  changes. 

6.3  mb  TEST  CALCULATIONS 

To  test  the  mb  algorithm  we  added  a  synthetic  seismogram  to 
seismic  noise,  scaled  to  different  amplitude  levels.  In  Figure  6.2 
we  showed  a  case  with  peak  synthetic/peak  noise  ■  100.  Examples 
will  be  shown  with  the  same  synthetic  and  the  ratio  reduced  to  3,  1 
and  0.5.  Two  different  noise  sections  were  used,  both  from  RKON 
recordings  that  were  in  the  AI  data  set.  The  timing  of  the  syn¬ 
thetic  within  the  noise  segment  was  also  varied. 

In  Figure  6.4  we  show  9  examples  that  are  variations  on  the 
example  in  Figure  6.2.  The  synthetic  seismogram  is  added  at  four 
different  times,  with  the  noise  scaled  to  different  levels.  In 
Figure  6.5  we  show  five  more  examples  which  are  constructed  the  same 
way  with  an  RKON  noise  sample  from  another  day. 


FILTER  CENTER  FREQUENCY  (Hz) 


Figure  6.3.  Comparison  of  the  Fourier  velocity  (amplitude 
times  frequency)  spectrum  to  the  MARS  log 
(A  '  f  )  spectrum  (dashed  line).  The  amplitude 
is  in  Sjj  units. 


Figure  6.5.  Test  of  the  m.  algorithm  with  noise  section  274110. 


Flltl-tMl  >11.  I  UP  *  »*M  I*1 

/  l  -.llil:  /  lt\U*:*4  JAII  1  /  llll  l  l;^': 

cmr  i  p  ii  /  i-  /!».•» *»»  /  i-t-i  »:  •*»/  i  >1/^1  m.  ; 


Figure  6.5.  (continued). 
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Figure  6.5.  (continued). 


fit  10 


(continued). 


The  "correct"  value  for  is  taken  from  the  peak  S/N  -  100 

cases.  For  the  62110  noise  sample  this  is  3.589  from  Figure  6.2. 
For  the  274110  noise  sample  it  is  3.586  from  the  first  plot  in 

Figure  6.5.  The  small  difference  is  due  to  differences  in  the  way 
the  synthetic  was  normalized  (note  the  amplitude  scale  on  the  peak 
S/N  -  100  synthetics). 

The  "correct"  arrival  time  is  assumed  to  be  the  tg  for  the 

Peak  S/N  -  100  cases.  It  differs  by  only  0.08  seconds  between  the 
two  noise  samples. 

The  important  data  from  these  test  calculations  are  summarized 
in  Table  6.1.  The  seismograms  are  identified  by  the  letter  given 

beside  the  m^  graph  on  each  plot.  We  also  list  the  S/N  ratio  at  1 

Hertz  and  the  The  m^  uncertainty  is  the  larger  of  the  DM3 

values,  which  is  always  the  negative  value.  It  depends  directly  on 
the  S/N  at  1  Hertz.  The  "m^  error"  is  simply  the  difference 

between  the  computed  m^  and  the  value  for  the  peak  S/N  *  100  case. 

The  magnitude  is  not  meaningful  if  it  is  not  computed  for  the 

right  signal.  Listed  in  the  table  is  the  "Arrival  Time  Error," 

which  is  the  difference  between  the  t„  for  the  peak  S/N  -  100 

9 

cases  and  the  tg  at  1  Hz  for  the  case  being  analyzed.  This  is 
also  indicated  on  the  plots.  The  lg  is  marked  with  an  asterisk. 
When  the  tg  at  1  Hertz  is  much  different  from  the  tg,  it  is 

marked  with  a  "+."  The  actual  arrival  time  of  the  signal  is  denoted 
by  tA,  as  explained  in  the  caption  to  Figure  6.4.  In  two  of  the 

low  S/N  cases,  I  and  N,  the  signal  was  entirely  missed.  In  those 
cases  an  arrow  marks  the  actual  arrival  time. 

There  are  six  cases  where  the  peak  S/N  -  1.  These  are  B,  C, 
0,  E,  K  and  L.  The  S/N  at  1  Hertz  is  2.0  to  2.6.  The  signal  can  be 
seen  visually  on  the  seismograms,  but  it  is  certainly  below  the 

threshold  where  one  woula  be  very  comfortable  with  time  domain  m^ 

measurements. 

For  each  of  the  peak  S/N  >  1  cases,  the  detection  algorithm 

accurately  located  the  signal.  The  maximum  arrival  time  error  for 
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TABLE  6.1 


SUMMARY  OF  m 5  TEST  CALCULATIONS 


Identifier 

Peak  Siqnal 
Peak  Noise 

S/N 
(1  Hz) 

A 

mb 

A 

Uncertainty 

E%or 

Arrival 

Time 

Error 

Noise  Section 

62110 

(Figure  6.4) 

(Figure  6.2) 

100 

217 

3.589 

0.001 

— 

— 

A 

3 

6.5 

3.583 

0.018 

-0.006 

+0.02 

B 

1 

2.2 

3.593 

0.055 

0.004 

0.03 

C 

1 

2.6 

3.667 

0.046 

0.078 

-0.09 

D 

1 

2.2 

3.586 

0.056 

-0.003 

0.48 

E 

1 

2.2 

3.581 

0.057 

-0.008 

-0.01 

F 

0.5 

1.5 

3.734 

0.082 

0.145 

0.09 

G 

0.5 

1.2 

3.622 

0.110 

0.033 

-0.76 

H 

0.5 

1.5 

3.708 

0.088 

0.119 

0.14 

I 

0.5 

1.4 

3.689 

Missed  Signal 

Noise  Section  274110  (Figure  6.5) 


J 

100 

219 

3.586 

0.001 

— 

- 

— 

K 

1 

2.0 

3.536 

0.111 

-0.05 

-0. 

.05 

L 

1 

2.2 

3.595 

0.095 

0.09 

-0. 

.15 

M 

0.5 

1.1 

3.583 

0.226 

-0.03 

-0. 

.05 

N 

0.5 

1.0 

3.548 

Missed  Signal 

the  1  Hertz  energy  is  0.48  seconds  for  0.  For  the  other  seven  the 
errors  are  0.15  seconds  or  smaller.  The  mean  error  is  0.04 

units,  with  three  of  the  six  having  errors  less  than  0.01  units. 

In  the  last  section  we  mentioned  some  algorithm  improvements 
that  should  be  added.  One  was  to  improve  the  method  for  selecting 
amplitude  peaks,  but  we  doubt  this  would  have  any  effect  on  the 

spectral  values  near  1  Hz  in  these  examples.  The  second  change  was 
to  fit  a  smooth  polynomial  to  the  spectral  values  between  0.8  and  1.2 
Hertz.  This  addition  would  probably  have  an  effect.  For  example, 
the  m^  for  K  would  increase,  reducing  the  error  for  that  case. 

There  are  also  six  cases  with  the  peak  S/N  =  0.5.  These  are 
cases  F,  G,  H,  I,  M,  N,  and  the  S/N  at  1  Hertz  is  1.0  to  1.5.  For 
these  the  signal  cannot  be  detected  by  eye.  Even  so,  the  1  Hz 

arrival  is  located  within  0.2  seconds  in  three  of  the  six  cases.  In 

one  case  the  error  is  0.76  seconds  and  there  are  two  cases  (I  ana  N) 
where  the  signal  is  missed  altogether. 

When  the  S/N  ratio  gets  very  low,  the  mb  errors  can  become 
significant,  as  should  be  expected.  For  example,  the  error  is  0.145 
units  for  Case  F.  We  must  be  careful  when  interpreting  this  error, 
however,  since  even  when  the  signal  is  missed  (Cases  I  and  N)  the 
computed  is  not  much  different  from  the  correct  value.  Of 
course,  it  shouldn't  be  since  the  noise  has  been  scalea  to  be  about 
the  same  size  as  the  signal. 

For  the  low  S/N  cases  the  t^  plot  is  Quite  discontinuous. 

The  spectral  amplitudes  are  being  selectea  from  different  portions 

of  the  signal.  For  example,  in  Case  F  the  1  Hertz  amplitudes  are 

associated  with  the  correct  arrival  time.  However,  the  1.6  Hertz 

amplitude  is  associated  with  an  arrival  time  28  seconds  later!  The 

t  -f  section  for  this  case  is  plotted  in  Figure  6.6.  There  are 
g  c 

several  apparent  undispersed  P  wave  arrivals,  indicated  by  a  series 
of  filter  peaks  aligned  at  the  same  time.  The  t  is  picking  one 
such  arrival,  ana  there  is  another  near  10530  on  the  t^  scale. 
The  actual  arrival  is  at  t,  and  the  algorithm  did  choose  the  1 

a 
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Hertz  amplitude  at  this  time.  However,  at  nearby  f  there  are 
spectral  amplitudes  associated  with  this  arrival  that  were  not 
chosen. 

The  point  of  the  discussion  in  the  last  paragraph  is  that  the 
t  -fc  plane  does  include  estimates  for  the  spectral  amplituae  of 
the  actual  signal  over  a  significant  bano  of  freauencies,  even  in 
low  S/N  cases  like  F.  Since  we  can  identify  the  correct  spectral 
peaks  by  eye,  the  automated  algorithm  to  pick  the  peaks  can  probably 
be  altered  to  duplicate  the  process.  This  would  improve  the  results 
for  very  low  S/N  arrivals.  But,  we  point  out  that  the  results  in 
Table  6.1  indicate  that  the  algorithm  is  already  giving  answers  that 
are  close  to  the  inherent  accuracy  threshold  imposed  by  the  noise, 
so  we  cannot  expect  too  much  improvement. 

6.4  tnb  FOR  RKON  RECORDINGS  OF  EURASIAN  EVENTS 

Using  the  algorithm  described  in  the  previous  section,  we 
computed  m^  for  RKON  recordings  of  eight  presumed  explosions  in 
the  Soviet  Union.  The  locations  of  these  events  are  shown  in  Figure 

6.7.  The  data  and  the  event  designations  are  from  the  AI  data  set 
(Savino,  et  aJL ,  1980). 

The  m^  determination  for  each  seismogram  is  shown  in  Figure 

6.8.  Note  that  the  signal /noise  at  1  Hertz  is  certainly  adequate 
for  accurately  computing  m^.  However,  there  would  be  less  noise 
contamination  if  the  m^  were  computed  at  a  higher  frequency  like 
1.2  Hertz.  There  is  no  reason  why  the  m^  should  not  be  computed 
at  some  frequency  other  than  1.0  Hertz,  as  long  as  it  is  done  consis¬ 
tently.  Since  the  process  is  automatic,  a  catalog  of  values  could 
be  maintained  for  future  examination. 

Similar  rri^  calculations  were  done  for  two  Tibetan  earth- 
duakes  recorded  at  RKON.  The  results  are  shown  in  Figure  6.9  where 
we  see  that  these  are  very  low  signal/noise  recordings,  in  fact,  the 
arrival  time  is  very  difficult  to  pick  by  eye.  Finally,  in  Figure 
6.10  we  show  the  m^  for  a  high  signal/noise  IASA  recording  of  a 
Kurile  arc  earthquake. 
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Figure  6.8.  mb  for  presumed  Soviet  explosions  recorded  at  RKON. 
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Figure  6.8.  (continued). 
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Figure  6.8.  (continued). 
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Figure  6.8.  (continued). 


Figure  6.8.  (continued). 
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Figure  6.8.  (continued). 
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Figure  6.8.  (continued). 


Figure  6.10.  m,  for  an  earthquake  recorded  at  LASA. 


In  Table  6.2  we  summarize  the  mb  data.  We  also  list  the 
signal/noise  at  1  Hertz  ana  the  statistical  uncertainty  in  m^. 
The  POE  mb  is  listed,  when  available.  Finally,  we  include  the 

mb  of  the  noise. 

The  noise  mb  is  useful  for  flagging  events  for  which  the 
calibration  appears  to  be  incorrect.  Events  274  and  276  are  in  this 
category.  Note  that  the  1  Hertz  noise  level  is  more  than  an  order 
of  magnitude  larger  than  it  is  for  the  other  RKON  recordings.  Also, 
for  276  the  m^  is  1.3  units  larger  than  the  POE  For  most  of 

the  other  events  the  mb  and  POE  mb  are  in  fairly  good  agree¬ 

ment.  The  main  exception  is  277,  and  the  calibration  may  also  be  in 
error  for  this  one. 

6.5  THE  M$  ALGORITHM 

The  m^  algorithm  is  a  natural  extension  of  the  MARS 

processing  which  leads  to  an  estimate  for  the  spectral  amplitude  of 
a  particular  signal  in  a  window  that  is  narrow  in  both  time  and 
frequency.  As  explained  in  Section  6.2,  the  "MARS  spectrum"  is 
constructed  by  selecting  an  envelope  peak  for  each  center 
frequency.  The  selection  is  done  by  finding  the  peaks  that  satisfy 
a  pre-selected  pattern.  For  mb  the  pattern  includes  the 
requirement  that  the  peaks  have  (approximately)  the  same  group 

arrival  time  (tg).  To  apply  the  same  algorithm  to  dispersed 
surface  waves,  it  is  again  necessary  to  specify  the  pattern  of  the 
signal  of  interest  in  a  display  like  that  in  Figure  6.1.  In  effect, 
the  program  must  be  given  some  estimate  for  the  group  velocity 

dispersion  of  the  signal  of  interest. 

If  the  tg  for  all  envelope  peaks  are  corrected  to  account 

for  a  specified  dispersion,  the  MARS  detection  algorithm  can  be  used 
just  the  same  as  it  is  for  body  waves.  A  corrected  display  of 

envelope  peaks  like  that  in  Figure  6.1  can  be  searched  for  an 
"undispersed"  arrival,  which  will  be  an  arrival  which  has  dispersion 
close  to  that  specified.  Once  the  arrival  is  detected,  M$  can  be 
computed  in  an  entirely  analogous  fashion  to  mb> 


126 


TABLE  6.2 


SUMMARY  OF  mb  FOR  ACTUAL  EVENT  RECORDINGS 


S/N 

/\ 

mw 

POE 

Noise 

Event 

Range 

(1  Hz) 

A 

mb 

Uncertainty  mb 

Presumed 

Explosions 

at  RKON 

20 

79.3 

92 

5.827 

0.002 

5.9 

3.86 

21 

71.0 

7 

5.114 

0.016 

4.8 

4.25 

22 

76.4 

6 

4.780 

0.025 

5.1 

4.03 

33 

53.8 

3 

4.622 

0.051 

5.0 

4.16 

81 

79.4 

58 

5.831 

0.003 

6.0 

4.07 

274 

76.5 

2 

5.466 

0.112 

not 

reported 

5.18 

276 

79.5 

7 

6.414 

0.027 

5.1 

5.59 

277 

79.4 

3 

4.326 

0.051 

5.0 

3.88 

Presumed 

Earthquakes 

at  RKON 

151 

93.6 

3.4 

4.922 

0.061 

5.0 

4.39 

159 

94.1 

1.6 

4.312 

0.105 

not 

reported 

4.10 

Presumed 

Earthquake 

at  LASA 

47 

61.1 

20 

5.116 

0.009 

4.9 

3.81 
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In  essence,  the  M$  algorithm  is  the  same  as  the  mb  algo¬ 
rithm,  except  that  we  must  first  map  the  tg-f  plane  into  a  new 
plane  in  which  signals  with  a  prespecifiea  dispersion  will  appear  to 
be  nearly  undispersed.  Masso,  et  al_.  (1979),  Appendix  C,  discuss 
this  mapping  and  suggest  two  ways  to  implement  it.  The  first  is  to 
select  some  reference  frequency  (f^)  and  compute 

«tg(f)  -  t*(f)  -  t*(fR)  ,  (6.1) 

where  the  t  *  are  the  group  arrival  times  for  the  prespecified 
dispersion.  The  t  -f  plane  is  then  mapped  into  a  t  1  — f  plane 
with  9 

t’(f)  «  tg(f)  -  «tg(f)  .  (6.2) 

The  second  method  suggested  by  Masso,  et  ah  (1979)  is  to 
cross-correlate  the  signal  with  a  reference  signal  and  then  to 
process  the  cross-correlation  time  series.  Such  a  "matchea 
filtering"  operation  to  improve  the  ability  to  estimate  surface  wave 
spectra  is  a  well-known  technique  which  is  incorporated  in  the 
COLLAPSE  program  at  Teledyne  Geotech  (Alexander  and  Lambert,  1971). 

In  most  cases  only  the  phase  of  the  reference  signal  can  be 
specified  with  any  confidence,  so  the  cross-correlation  reduces  to  a 
multiplication  by,  exp[-i$],  where  t  is  the  phase  of  the  reference 
signal . 

In  practice, 

4  -  -2wf  «t  (f)  (6.3) 

and  the  two  methods  only  differ  because  the  order  of  operations  is 
interchanged.  That  is,  in  the  first  method  the  signal  is  first 
narrow-band  filtered,  then  the  "matched  filter"  provided  by  the 
dispersion  of  the  reference  signal  is  applied.  In  the 
cross-correlation  method  the  "matched  filtering"  is  applied  before 
narrow-band  filtering. 


Once  the  t^-f  plane  has  been  mappea  to  a  t^’-f  plane,  the 

detection  algorithm  described  in  Section  6.2  can  be  applied  to 

A 

define  the  “MARS  spectrum"  of  a  detectea  signal.  The  is  then 

computed  from  the  smoothed,  noise-corrected  spectral  amplitude  at  a 
particular  period.  An  initial  choice  which  is  consistent  with 

current  practices  for  obtaining  a  time  aomain  M$  is  20  seconas, 
though  other  periods  may  turn  out  to  oe  better  when  large  amounts  of 
data  are  processed. 

While  the  procedure  outlined  appears  to  be  a  straightforwaro 
extension  of  the  mb  procedure,  there  are  several  important 
parameters  that  can  only  be  firmly  set  after  developing  some 

experience  with  actual  data.  For  example,  the  choice  of  Q,  which 
controls  the  width  of  the  narrow-band  filter,  involves  a  tradeoff 
between  the  desire  for  good  time  resolution  (for  detection)  and  good 
spectral  amplitude  resolution.  The  width  of  the  "acceptance  window" 
in  the  t  -f  plane  (Figure  6.1)  is  controlled  to  some  degree  by  an 
arbitrary  parameter  called  a  (Savino,  et  aU,  1980).  Small  a  impose 
the  requirement  that  the  signal  dispersion  closely  match  that  of  the 
reference  signal,  while  large  a  increase  the  possibility  that  the 
wrong  spectral  amplitudes  will  be  associated  with  the  signal. 

A 

Finally,  we  point  out  that  for  M^  we  only  need  to  detect  tne 
signal  in  a  narrow  period  oand,  initially  chosen  to  be  18  to  22 
seconds.  Thus,  large  differences  between  the  signal  aispersion  ana 
reference  dispersion  outside  this  band  only  matter  to  the  extent 
that  they  make  the  signal  hard  to  detect. 

How  do  we  propose  to  specify  the  reference  signal  dispersion? 
For  large  signal /noise  surface  wave  arrivals,  it  doesn’t  make  much 
difference,  the  correct  spectral  amplitudes  are  easily  identified. 
The  problem  is  with  small  signals.  The  reference  dispersion  must 
come  from  theoretical  models  or,  empirically,  from  larger  events  in 
the  same  source  region  as  the  event  of  interest. 

6.6  M$  TEST  CALCULATIONS 

A 

To  test  the  M$  calculation  and  optimize  the  selection  of  Q 
and  «,  we  followed  much  the  same  procedure  described  in  Section  6.3 
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for  m^.  That  is,  a  synthetic  seismogram  was  added  to  various 
levels  of  actual  seismic  noise.  In  this  case  the  noise  was  taken 
from  a  long  period  LASA  recording  in  the  AI  data  set.  Typical 
results  are  shown  in  Figure  6.11.  For  these  examples,  connection 
with  equations  (6.1  to  6.2)  of  the  previous  section  was  usea  to  map 
the  tg-f  plane  to  a  t^'-f  plane.  The  reference  signal 
dispersion  was  the  dispersion  used  to  compute  the  synthetic  itself. 

A 

In  Figure  6.11  we  first  show  the  Mg  calculation  for  a  signal 
that  is  essentially  just  the  synthetic.  The  noise  section  is 
included,  but  with  the  peak  noise  amplitude  scaled  to  be  1  percent 
of  the  peak  signal  amplitude  (peak  S/N  *  100).  Since  we  are 
subtracting  the  right  dispersion  for  this  signal,  the  t  is  close 
to  zero  over  the  entire  frequency  band.  The  deviation  from  zero  is 
due  to  the  inability  of  MARS  to  precisely  identify  the  group  arrival 
time  of  a  particular  frequency  component  and  some  errors  in  the 
“exact"  group  arrival  times  due  to  interpolation. 

The  other  two  examples  in  Figure  6.11  have  Peak  S/N  *  1  and 
have  been  added  at  different  times  in  the  noise  section.  The  data 
from  these  three  examples  are  listed  in  Table  6.3.  The  arrival  time 
error  is  based  on  the  tg  at  20  seconds,  which  is  listed  with  each 
example. 

The  examples  in  Figure  6.11  show  that  the  signal  is  accurately 

A 

detected  and  the  M$  values  accurately  represent  the  signal 
amplitude.  Note  that  the  M$  values  for  cases  B  and  C  are  nearly 
the  same,  while  the  peak  time  domain  amplitudes  are  much  different. 

Many  cases  like  those  in  Figure  6.11  were  run  with  peak  S/N  » 
1  and  2.  The  best  results  were  obtained  with  Q  =  250  f,  as  it  is 
for  the  examples  in  Figure  6.11,  and  Q  »  125f.  More  testing  needs 
to  be  done  to  choose  between  these  values. 

Tests  were  run  with  the  cross-correlation  method,  as  in 
Equation  6.3,  as  well  as  with  the  direct  mapping  of  the  tg-f 
plane.  Similar  results  were  obtained,  as  expected. 


r 


Figure  6.11.  (continued) 
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Figure  6.11  (continued) 


TABLE  6.3 
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SUMMARY  OF  Ms  TEST  CALCULATIONS 


Identifier 

Peak  Signal 
Peak  Noise 

S/N 

(20  sec) 

A 

A_ 

A 

Ms 

Uncertainty 

A 

Ms 

Error 

Arrival 

Time 

Error 

(sec) 

A 

100 

168 

2.169 

0.001 

— 

— 

B 

1 

3.3 

2.280 

.067 

.111 

-1.9 

C 

1 

4.4 

2.224 

.068 

.055 

3.8 
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6.7  CONCLUSIONS 

The  algorithm  is  essentially  ready  for  application  to  a 

large  data  base.  Some  minor  tinkering  with  the  algorithm,  as 
mentioned  at  the  end  of  Section  6.2,  might  improve  the  performance 
slightly. 

The  m^  data  will  not  be  compatible  with  standard  data. 
For  example,  8ache  (1979)  showed  that  the  m^-yield  curves  have 
different  slope  than  m^-yield  curves.  The  advantage  of  using 
is  that  it  Is  more  consistent  and  convenient  for  digital  data. 
However,  the  full  advantage  of  the  technique  requires  that  an  entire 
data  base  be  processed. 

A 

The  M$  algorithm  is  at  an  earlier  stage  of  development,  in 
part  because  we  do  not  have  the  large  base  of  experience  whicn  the 
VSC  discrimination  experiment  provided  for  short  period  body  waves. 

A 

However,  our  tests  so  far  are  promising  and  the  Ms  algoritnm  is 
clearly  ready  for  application  to  a  real  data  set. 

The  fi  will  never  be  as  automatic  as  mK  because  of  the 
s  b 

requirement  to  estimate  the  expected  dispersion  characteristics  of 
the  signal.  Of  course,  most  of  this  could  be  automated  as  part  of  a 
larger  signal  analysis  system. 
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